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ABSTRACT 


A  study  was  undertaken  to  examine  the  influence  of  manmade  changes  to  the 
drainage  network  on  flooding  from  urbanized  catchments.  The  literature  contains 
contradictory  hypotheses  regarding  the  influence  of  different  structural  changes.  Testing 
of  these  hypotheses  is  facilitated  by  the  use  of  a  physics  based  hydrologic  model.  A 
distributed  physically  based  hydrologic  model  Gridded  Surface  Subsurface  Hydrologic 
Analysis  (GSSHA)  was  combined  with  SUPERLINK,  a  storm  drainage  scheme.  Factors 
that  were  tested  include  storm  drainage  networks,  impervious  areas,  drainage  density, 
width  function,  and  initial  soil  moisture.  Results  indicate  that  the  subsurface  drainage 
network  is  important  for  moderate  rainfall  events,  but  largely  overwhelmed  for  extreme 
events.  Flood  magnitude  increases  due  to  modifications  to  the  channel  network  are 
dominant  compared  to  effects  from  impervious  area.  Initial  soil  moisture  likewise 
impacts  moderate  storms,  but  its  importance  diminishes  during  extreme  events. 


I.  INTRODUCTION 


There  is  no  argument  that  flood  magnitude  and  frequency  increase  as  urban 
development  spreads  throughout  a  watershed.  It  is  obvious  that  understanding  this  trend 
is  of  great  social  and  economic  importance.  But  what  causes  this  change  in  hydrology  is 
the  source  of  much  debate  and  numerous  studies.  The  hydrologic  processes  affected  by 
urbanization  are  primarily  infiltration  and  surface  runoff.  Comparisons  between  basins  of 
varying  land  use  and  drainage  systems  can  provide  some  insight  into  what  feature  plays 
the  dominant  role.  However,  it  is  impossible  to  eliminate  differences  in  scale, 
topography,  geometry,  and  geology  in  such  multi-basin  analyses.  A  numerical  model 
stands  out  as  the  best  method  to  systematically  isolate  watershed  properties  that  affect  its 
hydrologic  response. 

Distributed-parameter  models  are  gaining  a  foothold  in  the  world  of  hydrologic 
modeling,  as  these  methods  are  capable  of  describing  spatially  varied  land-surface 
modifications.  A  major  deficiency  of  most  these  models,  however,  is  their  inability  to 
explicitly  handle  storm  drainage  networks.  The  purpose  of  this  thesis  was  twofold:  to 
integrate  a  unique  storm  drainage  algorithm  to  an  existing  distributed  hydrologic  model, 
and  to  test  hypotheses  using  data  from  a  watershed  where  urbanizing  drainage  effects  are 
considerable.  This  study  consists  of  five  primary  hypotheses  relating  to:  the  storm 
drainage  network,  impervious  areas,  drainage  density,  the  distribution  of  drainage  within 
a  watershed,  and  pre-storm  soil  moisture. 

Changes  in  urban  runoff  volume  and  flood  peaks  have  historically  been  blamed 
on  increases  in  impervious  area.  This  theory  was  recently  challenged  by  a  study  in  and 
around  Charlotte,  North  Carolina  (Smith  et  al.  2002).  The  essence  of  their  conclusions 
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was  that  the  increase  in  storm  drainage  connectivity  and  hence  hydraulic  efficiency 
played  the  greatest  role  in  increasing  flood  magnitudes.  They  also  concluded  that 
antecedent  soil  moisture  and  drainage  density  distribution  within  the  watershed  play  an 
important  role  in  flooding. 

Storm  sewers  enhance  the  drainage  efficiency  of  a  watershed,  but  their  impact 
compared  to  other  modifications  is  unknown.  The  model  developed  here  provides  a 
unique  and  powerful  tool  for  assessing  hydrology  in  an  urban  catchment.  Storm  sewers 
are  often  overwhelmed  during  extreme  events,  so  it  is  possible  that  their  effect  diminishes 
as  the  precipitation  intensity  increases. 

Simulations  are  also  performed  with  different  combinations  of  impervious  area 
and  drainage  network  densities.  Analysis  of  these  results  identify  which  watershed 
feature  is  more  influential  in  increasing  flood  magnitudes.  Understanding  the  role  of 
impervious  coverage  versus  storm  sewers  could  be  useful  in  guiding  the  storm  water 
management  practices  of  urbanizing  areas,  as  well  as  solutions  to  flooding  in  urban 
environments. 

Drainage  density  represents  the  length  of  stream  channels  in  a  basin  relative  to 
total  drainage  area.  A  sensitivity  analysis  on  this  parameter  might  show  the  range  of 
drainage  densities  that  have  the  strongest  effect  on  flood  peaks.  This  dimensionless 
measurement  may  be  useful  to  predict  the  effect  of  development  on  other  watersheds. 
Similarly,  the  spatial  distribution  of  the  drainage  network  is  tested  for  its  influence  on 
basin  response.  This  distribution  is  measured  by  a  statistical  parameter  known  as  the 
width  function.  Relating  the  width  function  to  flood  peaks  likewise  might  provide 
generalized  conclusions  about  the  location  of  drainage  within  a  watershed. 
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Soil  moisture  preceeding  an  event  can  have  a  significant  impact  on  the  runoff 
generating  mechanisms.  Simulating  a  watershed  with  all  properties  being  identical 
except  initial  soil  saturation  points  to  the  importance  of  this  soil  parameter.  A  greater 
understanding  of  the  role  of  antecedent  soil  moisture  could  improve  flood  prediction 
accuracy.  A  coupled  surface-storm  sewer  model  has  the  potential  to  explore  each  of 
these  hydrologic  processes,  and  provide  conclusions  of  scientific  and  practical  impact. 
The  specific  objectives  of  this  thesis  are  to: 

1)  Implement  a  sophisticated  storm  network  model  and  couple  it  with  a  distributed 
hydrologic  model 

2)  Explore  the  importance  of  engineered  subsurface  storm  drainage  networks 

3)  Evaluate  the  role  of  impervious  areas 

4)  Compare  the  effect  of  varying  degrees  of  drainage  density 

5)  Determine  the  role  of  the  width  function  in  watershed  runoff 

6)  Evaluate  the  role  of  initial  soil  moisture  in  urban  basins. 
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II.  BACKGROUND 


Literature  Review 

Analysis  of  the  effects  of  urbanization  on  watershed  hydrology  is  not  new. 

Studies  dating  back  to  the  late  1800’s  have  been  concerned  with  flooding  in  settings  such 
as  urban  England.  It  is  therefore  critical  to  become  familiar  with  previous  attempts  to 
understand  the  nature  of  storm  drainage  systems  and  their  effect  on  flooding,  as  well  as 
general  changes  in  land  use.  The  availability  of  computing  power  has  revolutionized 
techniques  of  hydrologic  modeling,  and  thus  our  understanding  of  the  processes  has 
improved. 

There  is  no  question  that  land  surface  modifications  related  to  urban  development 
have  increased  the  magnitude  and  frequency  of  flooding  around  the  globe.  Each  of  the 
research  efforts  included  in  the  following  review  cites  the  importance  of  understanding 
the  nature  of  anthropogenic  effects  on  a  catchment’s  drainage  characteristics.  These 
observed  changes  in  stream  flow  have  been  modeled  with  empirical  techniques  largely 
dependent  on  assumptions  of  the  controlling  factors  of  urban  runoff  Indeed,  the 
commonly  accepted  principal  variable  for  much  of  the  20*  century  has  been  the  extent  of 
impervious  area.  Eeopold  (1968)  presented  an  empirical  analysis  of  impervious  coverage 
that  became  a  benchmark  of  urban  watershed  theory.  Essentially,  this  research  showed 
that  increases  in  impervious  area  were  accompanied  by  elevated  flood  peaks  and  a 
decrease  in  inter-storm  flows.  Although  Eeopold  acknowledged  the  concept  that  storm 
sewers  decrease  the  lag  time,  the  drainage  network  was  given  a  secondary  role  to  land 
use.  Eeopold  compiled  a  table  based  on  contemporary  research  of  Carter  (1961), 
Anderson  (1968)  and  others.  By  providing  two  values,  percentage  impervious  and 
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percentage  sewered,  one  could  estimate  the  increase  in  flood  peaks.  Based  upon  this 
theory,  simplifying  assumptions  have  been  made  in  urban  hydrologic  modeling. 

Primarily,  it  has  been  assumed  that  the  effect  of  urbanization  can  be  described  by  the 
fraction  of  impervious  area.  Secondly,  it  is  commonly  thought  that  large-scale  events 
will  overwhelm  a  storm  drainage  system,  and  thus  their  total  contribution  to  runoff  is 
small. 

Climate 

Many  studies  have  focused  on  other  reasons  for  increasing  flood  magnitudes  and 
frequencies.  Some  have  pointed  to  changing  climate,  including  work  by  Reynard  et  al. 
(2001).  These  studies  mainly  pertain  to  large  catchments  on  the  order  of  10,000  sq  km, 
and  are  far  too  coarse  to  describe  the  small-scale  effects  of  urbanization.  Howe  and 
White  (2000)  refuted  the  claim  that  global  warming  is  to  blame  for  increased  urban 
flooding.  From  a  political  policy  standpoint,  they  criticized  the  practice  of  straightening 
rivers  and  expanding  local  drainage  systems.  This  common  solution  to  localized  flooding 
concerns  simply  moves  the  problem  downstream,  as  both  time  to  peak  and  attenuation  are 
decreased.  They  note  that  there  must  be  a  comprehensive  plan  that  incorporates  all  scales 
of  drainage  in  order  to  reverse  the  current  trend  of  flooding.  Their  suggestions  include 
constructing  wetlands  to  recreate  the  original  flow  attenuation. 

Land  Use 

Land  use  changes  have  been  the  focus  of  many  flood  related  studies.  Beighley 
and  Moglen  (2003)  attempted  to  compare  flood  peaks  from  pre-urbanized  flow  records  to 
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current  urbanized  flow  data.  For  example,  they  postulate  that  a  100  year  flood  in  the 
1950’s  would  be  a  75  yr  flood  in  the  1990’s.  Their  method  of  correcting  flows  to 
represent  statistieally  similar  floods  was  applied  to  current  watershed  conditions.  They 
focused  on  three  factors:  increases  in  urban  land  use,  decreases  in  forested  land,  and  the 
percent  overall  urbanization  classified  as  “high  density  development”.  It  was  found  that 
none  of  these  land  use  factors  were  effeetive  for  estimating  the  adjusted  flood  flows. 
Surprisingly,  a  spatial  parameter  describing  the  loeation  in  the  basin  in  which 
development  oeeurred  could  best  prediet  the  adjusted  flows.  Development  oceurring 
farthest  from  the  basin  outlet  had  the  greatest  effect  on  peak  flows,  and  thus  they  noted 
the  need  to  investigate  ehanges  in  channel  roughness.  The  use  of  a  detailed  hydrologic 
model  could  explore  both  the  land  use  question  and  spatial  observation.  Crooks  and 
Davies  (2001)  presented  additional  researeh  on  the  issue  of  land  use.  Using  30  years  of 
data  of  land  use  change  and  its  flooding  effeets  in  the  Thames,  England  catchment,  they 
found  that  eoverage  had  an  insignificant  role  compared  to  rainfall  at  a  large  scale.  These 
findings  point  to  the  fact  that  another  proeess  dominates  flood  dynamies  in  an  urbanizing 
watershed. 

Channel  &  Drainage  Network  Morphology 

Researeh  into  the  role  of  drainage  networks  in  flooding  is  not  new.  Some  of  the 
earliest  work  was  done  by  Anderson  (1970)  on  the  effects  of  urban  development  on 
floods  in  northern  Virginia.  He  gathered  rainfall  and  runoff  data  on  8 1  similar 
watersheds,  developed  relationships  between  size,  length,  slope,  pereent  impervious,  and 
type  of  drainage  system.  Using  a  comparative  approach  between  basins  with  one  or  more 
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similar  characteristics,  Anderson  predieted  the  effeet  of  ehanging  one  variable.  Drainage 
system  improvements  were  shown  to  reduce  the  lag  time  of  a  flood  by  up  to  12%.  The 
same  lag  time  was  predieted  for  a  sewered  system  with  impervious  areas  versus  pervious 
surfaees,  but  with  a  higher  magnitude  flood.  Regardless  of  land  surfaee  ehanges,  the 
effeet  of  sewer  installation  was  to  inerease  the  flood  peak  magnitude  by  a  faetor  of  2  or  3 
for  smaller  events.  While  a  eompletely  impervious  surfaee  eould  inerease  the  flood 
magnitude  by  2.5  times,  the  effect  of  this  land  surface  eondition  deereased  as  the  storm 
total  rainfall  grew.  In  faet,  he  predieted  that  the  land  surfaee  modifieations  in  some 
catchments  were  negligible  for  a  100-year  flood,  as  the  entire  watershed  would  behave  as 
impervious.  Though  these  results  were  based  on  estimates  and  empirieal  relationships, 
the  eonelusions  certainly  opened  the  way  for  new  drainage  network  studies. 

In  more  reeent  ehannel  researeh,  Wolff  and  Burges  (1994)  used  the  model 
DAMBRK  to  explore  the  effeets  of  river  channel  properties  on  downstream  flood 
frequeney.  They  eoneluded  that  an  increase  in  storage  in  the  channels  deereases  the  flood 
peaks.  An  interesting  possibility  elueidated  by  Wolff  and  Burges  (1994)  is  that 
alternatively,  eovered  channels  (culverts)  deerease  the  storage  eapaeity,  foreing  exeess 
flow  to  the  overland  plain,  elevating  loeal  flooding. 

Despite  the  focus  of  this  study  on  urban  and  developing  watersheds,  researeh  on 
an  undeveloped  watershed  by  Troeh  et  al.  (1994)  introdueed  relevant  findings  related  to 
flow  veloeities.  Using  empirieal  and  model-based  analysis  from  a  small  eatehment  in  the 
Appalaehians,  they  sought  to  draw  a  eonneetion  between  peak  diseharge  times  and  flow 
veloeities.  For  lower  overland  flow  veloeities  (due  to  high  roughness  values),  the 
ehannel  veloeity  did  not  infiuenee  flood  peak  timing.  However,  as  overland  roughness 
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decreased,  the  channel  properties  began  to  control  lag  times.  Below  a  critical  overland 
roughness  value,  time  to  peak  was  completely  governed  by  channel  roughness.  Although 
this  critical  overland  roughness  threshold  was  below  a  naturally  occurring  value,  the 
study  suggests  that  urban  land  covers  with  low  roughness  (i.e.  rooftops  and  parking  lots) 
transfer  control  of  lag  time  from  the  overland  plain  to  the  channel  properties. 

Related  modeling 

Much  of  this  review  has  been  of  conceptual  and  regional  scale  studies.  Since  the 
focus  of  this  research  is  on  small  scale  distributed  hydrologic  modeling,  a  review  of 
similar  modeling  efforts  is  relevant.  A  comparable  approach  to  storm  drainage  modeling 
was  found  in  the  literature  by  Hsu  et  al.  (2000).  A  heavy  typhoon  was  simulated  over  a 
catchment  subject  to  flooding  in  Taipei,  Taiwan.  Much  of  the  city  is  protected  by  levees 
and  thus  dependant  on  pump  stations  to  remove  storm  water.  The  focus  of  their  study 
was  on  areas  of  inundation,  not  considering  the  effect  of  the  drainage  network  on  lag  time 
and  flood  magnitudes.  Using  the  storm  sewer  module  of  SWMM  and  a  two-dimensional 
diffusive  wave  overland  flow  routine,  the  authors  showed  that  the  storm  drains  increased 
the  areas  of  ponded  water.  Linkage  of  the  two  models  was  not  well  explained,  but  it 
appeared  that  the  models  were  not  run  simultaneously.  In  a  situation  where  flooded 
network  discharge  back  onto  the  overland  plain  is  a  consideration,  it  is  critical  that  the 
model  components  simulate  the  hydrologic  processes  in  parallel. 

Additional  modeling  to  assess  the  causes  of  increasing  flood  peaks  near  Charlotte, 
North  Carolina  Turner-Gillespie  et  al.  (2003)  focused  on  attenuating  reaches.  In  this 
study,  floodplain  roughness  was  tested  for  its  control  over  flood  peaks  and  lag  timing. 
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Although  decreasing  roughness  in  overland  flow  associated  with  urbanization  did 
produce  higher  peaks  and  quicker  lag  times,  the  relative  effect  was  small.  However,  it 
was  found  that  watershed  scale  flood  peaks  were  very  sensitive  to  the  tributary  response 
time.  This  same  conclusion  could  likely  be  drawn  for  storm  drainage  networks,  which 
generally  replace  natural  tributaries. 

Width  function 

An  investigation  into  the  nature  of  urbanizing  stream  networks  by  Graf  (1977) 
showed  a  drastic  increase  in  the  number,  length,  and  density  of  man-made  channels  in  an 
Iowa  catchment  with  extensive  historical  data.  Though  he  did  not  specifically  address 
subterranean  conduits,  the  research  on  artificially  altered  streams  is  still  applicable. 
Statistical  analysis  with  power  law  equations  showed  a  high  sensitivity  of  time  to  peak 
and  kurtosis  (peakedness)  for  channel  modifications. 

The  density  of  links  in  a  drainage  network  mentioned  by  Graf  (1977)  can  be 
described  by  a  statistic  now  known  as  the  width  function.  Generally,  the  width  function 
is  defined  as  a  plot  of  the  number  of  channel  segments  at  a  specified  distance  from  the 
basin  outlet  (Rodriguez -Iturbe  and  Rinaldo,  1997).  Attempts  have  been  made  to  relate 
the  width  function  to  peak  discharges,  including  work  by  Veitzer  and  Gupta  (2001). 
Focusing  strictly  on  network  properties,  they  determined  that  the  width  function  alone 
does  not  provide  enough  to  formulate  a  substantial  connection  between  drainage  structure 
and  flow  peaks.  It  is  clear  that  additional  information  is  necessary  for  the  width  function 
to  be  useful.  A  slightly  different  approach  defines  the  width  function  as  the  drainage  area 
at  a  flow  distance  from  the  outlet  (Richards-Pecou,  2002).  Although  their  particular 
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study  assumed  a  constant  flow  velocity  in  the  channels,  this  work  showed  increasing 
peaks  as  the  width  function  was  elevated  over  natural  conditions. 

An  excellent  example  using  the  width  function  concept  was  performed  by  Smith 
et  al.  (2002)  who  studied  the  heavily  urbanized  Charlotte,  North  Carolina  basin.  Five  of 
the  largest  flood  peaks  in  a  74  year  record  had  occurred  within  the  last  6  years  of  record 
on  the  Little  Sugar  Creek  catchment.  In  an  attempt  to  relate  the  extensive  urban  and 
suburban  growth  since  1960  to  the  increased  flooding,  their  research  examines 
hydrometeorology,  soil  moisture,  impervious  area,  and  drainage  network  modification. 
The  methodology  concentrates  on  diagnostic  testing  of  each  large  event,  as  the  target  was 
regional  hydrology.  However,  their  findings  provide  hypotheses  for  our  study.  A 
detailed  physical  model  can  be  used  to  investigate  the  relative  effect  of  these  hydrologic 
and  hydraulic  elements. 

Changes  in  landuse  have  traditionally  been  blamed  for  flooding,  but  surprisingly 
did  not  play  a  major  role  in  overall  water  balance  for  the  five  extreme  events.  The  most 
important  landuse  type,  impervious  area,  did  correlate  to  flood  peak  timing  and 
magnitude.  However,  the  extreme  events  of  interest  consist  of  rainfall  rates  far  higher 
than  saturated  hydraulic  conductivity  values.  Thus  infiltration  excess  is  the  dominant 
runoff  mechanism  with  or  without  the  impervious  coverage.  Antecedent  soil  moisture 
was  shown  to  play  an  important  role  for  some  events,  as  dry  soil  conditions  consumed  an 
appreciable  amount  of  water  through  infiltration. 

Of  greatest  interest  is  their  conclusion  that  expansion  of  the  drainage  network 
played  a  central  role  in  the  rising  trend  of  flood  peaks.  A  distinct  alteration  in  width 
function  could  be  seen  in  the  lower  section  of  Little  Sugar  Creek,  having  the  effect  of 
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decreasing  lag  time  and  increasing  flood  magnitudes.  Attenuating  reaches  were  shown  to 
have  significant  impact  on  reducing  peak  discharge,  and  conversely  channelizing  these 
reaches  would  have  the  opposite  effect. 
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storm  Drainage  Model  Selection 


GSSHA 

The  need  to  simulate  surface  water  flows  in  watersheds  with  diverse  runoff 
production  mechanisms  has  prompted  the  development  of  a  physically  based  hydrologic 
model,  called  the  Gridded  Surface/Subsurface  Hydrologic  Analysis  (GSSHA)  model 
(Downer  and  Ogden,  2004).  GSSHA  simulates  stream  flow  generated  by  both  infiltration 
excess  and  saturation  excess  mechanisms.  The  model  employs  mass-conserving 
solutions  of  partial  differential  equations  and  closely  links  the  hydrologic  components  to 
assure  an  overall  mass  balance.  It  is  not,  however,  capable  of  modeling  storm  sewers.  A 
sub-model  within  the  GSSHA  framework  to  directly  handle  subterranean  sewers  and 
grate  inlets  would  allow  modeling  of  complex  urban  catchments.  The  concept  of 
numerical  storm  water  modeling  is  not  new,  and  thus  the  selection  of  a  modeling 
approach  required  review  of  existing  methods.  This  review  includes  commonly  available 
models  as  well  as  some  lesser-known  methods,  and  a  short  description  of  each. 

uses  Full  Equations  model  (FEQ) 

The  FEQ  model  (USGS,  1997)  has  been  shown  to  accurately  model  free  surface 
flow  using  a  four  point  implicit  Preissmann  scheme.  For  conduits  flowing  full,  the  free 
surface  assumption  disappears.  This  problem  is  commonly  avoided  by  the 
implementation  of  a  “Preissmann  slot”,  a  narrow  slot  extending  from  the  top  of  the  pipe 
in  which  flow  is  permitted.  The  weakness  of  this  model  is  that  the  slot  is  assumed  never 
to  be  overtopped,  thus  preventing  surcharged  manholes  or  flow  back  into  the  overland 
plain.  Additionally,  Meselhe  and  Holly  (1997)  demonstrated  that  for  Froude  numbers 
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near  1.0,  the  Preissmann  4  point  scheme  is  unstable.  This  creates  problems  where  bed 
slopes  change  rapidly  and  other  trans-critical  flow  situations  as  well  as  free  jet  outlet 
conditions.  Given  these  limitations,  the  FEQ  model  was  not  found  suitable  for  general 
storm  sewer  modeling. 

U.S.  National  Weather  Service  DWOPER 

The  Dynamic  Wave  OPERational  (DWOPER)  model,  developed  by  Eread 
(1976),  solves  the  conservation  of  mass  and  conservation  of  energy  equations.  This  one 
dimensional,  unsteady  flow  model  allows  for  wave  propagation  in  the  upstream  and 
downstream  direction  to  account  for  situations  of  hurricane  surge  and  complex  backwater 
effects.  DWOPER  is  limited  in  cross  section  interpolation  and  its  inability  to  model 
supercritical  flow  (NWS,  2003). 

SWMM 

Development  of  the  U.S.  Environmental  Protection  Agency’s  Storm  Water 
Management  Model  (SWMM)  began  in  the  late  60’s,  and  has  since  undergone  many 
improvements.  The  model  contains  four  distinct  modules.  Runoff,  Transport,  Extran,  and 
Storage/Treatment.  The  Extran  module  was  of  interest  to  this  review,  as  it  allows  for 
pipe  flow  with  backwater,  surcharge,  pressurized  flow,  and  looped  networks.  Solving  the 
de  St-Venant  equations  explicitly  is  a  noted  weakness,  as  it  tends  to  become  unstable  at 
long  time  steps  and  is  computationally  inefficient.  Because  the  existing  code  is  a  product 
of  30  years  of  development,  it  would  be  difficult  to  merge  a  specific  component  to  suit 
the  needs  of  this  research  (USEPA,  2002).  Although  the  planned  overhaul  of  SWMM 
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would  make  it  more  adaptable  to  other  applieations,  there  are  eomponents  of  the  model 
sueh  as  the  land  surface  parameterization  scheme  that  cannot  match  the  sophistication  of 


GSSHA. 


Danish  Hydraulics  Institute  MOUSE 

The  MOUSE  pipeflow  model  accounts  for  all  of  the  critical  elements  of  a  storm 
drainage  model  including  flow  reversal,  backwater  effects,  pressurized  flow,  surcharged 
manholes,  and  storage  basins.  It  solves  the  full-dynamic  form  of  the  de  St.  Venant 
equations  by  an  implicit,  finite  difference  formulation  with  an  adaptive  time  step. 
However,  a  thorough  analysis  of  the  techniques  employed  in  MOUSE  is  impossible.  The 
source  code  is  not  available  since  MOUSE  is  proprietary  software,  making  re¬ 
development  too  complex. 

SUPERLINK  Scheme 

The  SUPERLINK  model,  developed  by  Ji  (1998),  is  a  general  hydrodynamic 
model  for  sewer/channel  networks.  This  method  solves  the  full  dynamic  de  St.  Venant 
equations  in  one  dimension  and  employs  the  Preissman  Slot  to  extend  the  open  channel 
flow  assumptions  to  closed  conduits  flowing  full  and  surcharged.  Unlike  many 
applications  of  the  Preissman  slot,  SUPERLINK  does  not  consider  the  area  of  the  narrow 
slot  in  the  flow  calculations  or  wetted  perimeter,  thus  reducing  the  error  associated  with 
this  concept.  Another  significant  benefit  of  this  model  is  its  staggered  grid  implicit 
solution  to  the  full  equations  of  motion,  enhancing  stability  and  computational  speed. 
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Of  these  sehemes,  SUPERLINK  was  judged  the  most  eapable  with  a  very  detailed 
description  of  the  formulation  allowing  implementation.  Although  it  has  not  had  the 
widespread  use  and  acceptance  as  many  of  the  other  models,  Ji  (1998)  tested  this  scheme 
on  a  complex  data  set  from  the  city  of  Winnipeg,  Manitoba,  Canada.  Winnipeg,  located 
on  the  banks  of  the  Red  River,  is  a  very  low  gradient  watershed  (subject  to  backwater 
effects  and  surcharging)  and  contains  multiple  looped  and  branched  pipes.  Ji  compared 
SUPERLINK  output  to  both  SWMM  Extran  and  physical  observations,  with  favorable 
results  against  both.  SUPERLINK  was  stable  at  a  400  second  time  step,  compared  to  a  7 
second  Extran  time  step,  and  generated  mass  conservation  errors  of  only  0.32%  over  a  5- 
hour  simulation  period.  Based  on  the  comprehensive  formulation  and  satisfactory 
simulation  results,  SUPERLINK  was  chosen  as  the  pipe  network  model  to  couple  with 
GSSHA. 
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III.  PIPE  NETWORK  DEVELOPMENT  &  TESTING 


SUPERLINK  Review 

Network  Nomenclature 

To  understand  the  following  discussion,  it  is  important  to  review  the  network 
definitions  and  rules.  Superlinks  are  series  of  links  connecting  junctions,  and  must  have  a 
junction  on  either  end.  A  junction  is  defined  as  a  point  where  two  or  more  superlinks 
meet,  or  the  unconnected  end  of  a  superlink  (such  as  intake/discharge  point  of  network). 

A  link  is  a  segment  of  a  superlink  connecting  two  nodes,  and  a  node  is  a  computational 
point  in  a  superlink.  Nodes  exist  adjacent  to  both  upstream  and  do wnsteam  junctions  of  a 
superlink,  thus  the  number  of  nodes  is  equal  to  one  plus  the  number  of  links. 

Each  junction  is  assigned  a  number,  1  through  number  of  junctions,  generally 
ordered  from  top  of  catchment  down.  Each  superlink  is  assigned  a  number,  1  through 
number  of  superlinks,  in  general  order  from  top  of  catchment  down.  Within  a  superlink, 
node  numbering  starts  at  the  first  upstream  node  (adjacent  to  upstream  junction)  and 
continues  through  number  of  links  plus  1 .  A  link  has  the  same  index  as  its  upstream 
node,  1  through  number  of  links.  The  nature  of  drainage  networks  often  includes  looped 
systems,  and  thus  the  numbering  scheme  must  be  capable  of  representing  this.  Although 
the  algorithm  is  not  limited  by  the  order  in  which  components  are  labeled,  computational 
efficiency  is  gained  through  organized  numbering.  Eor  model  stability  and  accuracy,  a 
minimum  number  of  3  links  (4  nodes)  are  required  for  each  superlink. 
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Superlink  1 


Junction  1 


Assumed  Flow  Direction 


Junction  2 


Figure  1  SUPERLINK  junction,  link,  and  node  nomenclature 
Inflow  is  allowed  at  junctions  and  nodes  via  two  primary  structures.  The  first  is  a 
culvert,  which  captures  a  natural  stream  channel,  and  only  is  possible  at  a  junction.  The 
second  is  any  type  of  grate/curb  opening  in  a  roadway,  and  is  possible  at  either  a  junction 
or  a  node.  Discharge  can  occur  from  either  a  flooded  manhole,  drop  inlet,  or  an  outlet 
pipe  (node  or  junction),  and  junctions  may  discharge  directly  into  a  channel.  The 
parameters  for  these  inlets  and  outlets  are  contained  in  the  node  and  junction  cards 
defined  in  the  file  format. 

The  use  of  both  nodes  and  links  may  seem  redundant,  but  in  fact  is  quite  integral 
to  the  “staggered  grid”  technique  employed  in  SUPERLINK.  There  are  two  variables 
that  must  be  solved  throughout  the  system  at  each  time  step;  depth  and  flow.  In  this 
approach,  depth  and  flow  are  considered  at  two  distinct  locations;  depth  is  computed  at 
nodes,  and  flow  in  links.  Calculating  these  two  variables  at  the  same  point  is  a  source  of 
instability  because  they  are  not  independent.  SUPERLINK  is  capable  of  modeling  long 
time  steps  due  to  this  staggered  feature  coupled  with  the  implicit  solution  to  be  discussed. 


Modeling  Theory 

The  basic  concepts  of  the  SUPERLINK  scheme  have  been  briefly  outlined.  Here, 
the  theory  behind  the  algorithm  will  be  explored  further.  The  central  equations  solved  in 
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this  model  are  the  conservation  of  mass  (1)  and  the  de  St.  Venant  equation  of  motion  (2). 
This  pair  of  nonlinear  partial  differential  equations  take  the  form  of: 


where  A  =  flow  cross-sectional  area,  Q  =  discharge,  h  =  depth,  u  =  velocity.  So  =  bed 
slope  of  conduit,  Sf=  friction  head  loss  slope,  Sl  =  local  head  loss  slope,  qo  =  lateral  flow 
to  conduit,  g  =  gravitational  constant,  x  =  distance,  and  t  =  time. 

The  two  central  equations  (1  &  2)  are  applied  on  sections  of  a  conduit  segmented 
by  computational  nodes.  Conservation  of  mass  is  represented  by  Equation  1,  and  is 
applied  across  a  node.  The  staggered  grid  approach  requires  the  conservation  of 
momentum  Equation  2  to  be  applied  on  a  different  control  volume.  The  layout  of  these 
volumes  is  shown  by  the  following  graphic. 


Node  1 


Node  2  Link  2  Node  3 


Link  1 


Eigure  2  SUPERLINK  staggered  grid  computational  scheme 


The  control  volume  shown  in  short  dashes  illustrates  the  continuity  equation  for  node  2, 
while  the  long  dashed  envelope  indicates  the  momentum  equation  for  link  2. 

The  St.  Venant  equations  of  motion  only  apply  to  free  surface  flow.  This  fact 
limits  their  application  to  pressurized  flow,  but  switching  to  a  closed  conduit  equation  is 
often  an  unstable  transition.  As  mentioned  in  the  model  review  section,  a  common 
solution  is  to  employ  the  “Priessmann  slof  ’  to  extend  the  free  surface  equations  to 
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conduits  flowing  full.  Again,  this  slot  area  is  not  used  for  flow  ealeulations,  but  merely 
to  pressurize  the  conduit  still  being  modeled  by  open  channel  flow  equations. 


Figure  3  Cross  section  view  of  pipe  with  Priessmann  Slot 


Linearized  Equations 

To  solve  the  partial  differential  equations,  they  must  be  discretized  over  their 
respective  control  volumes.  Thus,  unsteady  terms  such  as  flow  rate  and  depth  beeome 
time  dependent  variables.  The  diseretized  continuity  equation  with  indiees  referring  to 
Fig.  2  becomes 
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The  momentum  equation  takes  a  similar  form 


(e;*"  _  g,'  I  ^ )  +  +  gA^  +  S,,  K  =  +  gA^  -  A, 


At 


(4) 


The  only  new  term  in  these  equations  is  B,  or  the  top  width  of  flow  area. 
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Boundary  Conditions 

As  with  any  modeling  problem,  a  set  of  boundary  conditions  must  be  applied  to 
the  extents  of  the  network.  With  regards  to  the  SUPERLINK  model,  these  boundaries  are 
located  at  the  ends  of  each  superlink,  or  junctions.  The  first  component  of  the  junction 
boundary  is  the  water  surface  elevation  (head).  Junction  heads  may  be  known  or 
unknown,  as  determined  by  the  actual  network  configuration.  A  known  junction  head 
may  be  controlled  by  something  external  to  the  model,  such  as  a  reservoir  at  the  network 
outlet.  This  feature  would  create  backwater  pressure  propagating  upstream,  thus 
affecting  flow  upstream.  Unknown  junction  heads  occur  at  internal  connections  of  two  or 
more  superlinks.  Junctions  representing  an  intake  structure  at  the  start  of  a  superlink 
could  also  have  an  unknown  head. 

Flow  into  and  out  of  these  junctions,  whether  of  known  or  unknown  head,  is 
governed  by  end  condition  boundary  equations.  Inlet  entrance  geometry  governs 
pipeflow  in  steep  channels,  and  exit  properties  can  control  in  low  gradient  conditions. 

The  end  equations  use  the  head  in  the  junction  as  well  as  geometric  variables  to  produce  a 
set  of  coefficients  for  each  inlet  and  outlet.  The  inlet  and  outlet  coefficients  by  Ji  (1998) 
were  found  to  be  unstable  in  certain  situations  and  were  reformulated  as  discussed  in  the 
model  development  section. 

Solution  Technique 

The  implicit  scheme  is  defined  by  a  simultaneous  solution  to  all  unknowns  in  the 
system  at  each  time  step.  Instead  of  computing  the  head  at  every  internal  point  (junctions 
and  nodes)  as  the  model  steps  through  time,  only  unknown  junctions  are  part  of  the 
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solution  matrix.  The  reduction  in  the  matrix  size  and  thus  computational  demand  is 
substantial.  But  the  elegance  of  this  routine  is  the  way  in  which  the  unknown  internal 
node  depth  and  flow  are  incorporated  into  the  junction  matrix  solution.  Through  a  series 
of  recurrence  relations,  the  momentum  and  continuity  equations  are  propagated 
throughout  each  superlink  from  one  node  to  the  next.  This  is  done  in  both  the  forward 
and  reverse  directions  in  order  to  capture  both  positive  and  negative  flow.  The  resulting 
coefficients  become  part  of  a  relatively  complex  equation  relating  junction  heads, 
superlink  end  conditions,  internal  node  depth,  internal  pipe  flow  rate,  and  current 
timestep  network  inputs. 

The  matrix  used  to  solve  these  equations  takes  the  form  as  shown,  and  typically 
could  exceed  200  x  200  in  size. 
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In  this  matrix,  three  main  components  are  represented.  First,  the  square  matrix  set  is 
clearly  of  a  sparse  nature,  but  randomly  so.  The  diagonal  elements  denoted  by  F  and  off- 
diagonal  elements  O  and  \|/  all  are  summations  of  superlink  coefficients  (equations  30-33, 
Ji,  1998).  The  particular  coefficients  summed  as  well  as  their  location  in  the  matrix  are 
dependent  on  the  connectivity  of  the  network.  For  example,  a  junction  with  only  one 
superlink  to  or  from  would  only  see  one  off  diagonal  element,  whereas  a  complex 
junction  with  multiple  connections  would  have  numerous  off-diagonals.  The  second 
variable,  H,  is  a  one-dimensional  array  of  heads  at  the  future  time  step  for  each  junction 
in  the  network.  Third  is  the  “right  hand  side”  vector  of  the  matrix  equation.  G  is 
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composed  of  values  from  the  eurrent  time  step  ineluding  junetion  head  and  boundary 
conditions,  as  well  as  junction  inflows  for  the  future  time  step. 

With  the  matrix  eonstrueted,  any  number  of  solution  teehniques  ean  be  applied. 
The  order  of  the  left  hand  square  matrix  is  equivalent  to  the  number  of  junctions  in  the 
system.  If  the  order  beeomes  large  (over  500)  it  would  be  prudent  to  consider  efficient 
solving  routines  designed  to  utilize  the  sparse  nature  of  the  matrix.  The  location  of  non¬ 
zero  elements  is  eritieal  to  redueing  eomputation  time.  Although  the  majority  of  non¬ 
zero  elements  in  this  case  will  be  diagonally  biased,  there  are  cases  where  an  outlier 
would  render  some  solution  teehniques  useless.  Thus  the  generalized  LU  decomposition 
solution  method  was  chosen,  which  provides  a  full  solution  of  a  sparse  matrix. 
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Model  Development 


SUPERLINK  Algorithm 

The  actual  coding  of  the  SUPERLINK  scheme  initially  appeared  to  be  a  trivial 
task.  Ji’s  (1998)  publication  was  carefully  laid  out  to  provide  a  step-by-step  review  of  the 
model’s  formulation.  The  task  of  coding  was  undertaken  in  the  computer  language  C.  It 
soon  became  apparent  that  there  were  unforeseen  difficulties  in  reproducing  a  functional 
code.  The  very  nature  of  this  implicit  staggered  method  can  be  likened  to  an  automatic 
transmission,  contrasted  to  the  “manual  shift”  finite  volume  concept.  Although  this  may 
seem  humorous,  nothing  was  obvious  about  the  numerics  of  SUPERLINK.  It  was 
extremely  difficult  to  trace  bugs  in  the  developing  code  because  numerical  instabilities 
seemed  to  propagate  from  nowhere.  Ji’s  (1998)  algorithm  was  reproduced  by  following 
these  general  steps; 

1)  Calculate  momentum  and  continuity  coefficients  for  a  superlink 

2)  Calculate  forward  recurrence  relations  for  each  node  and  link  within  the  superlink 

3)  Calculate  reverse  recurrence  relations  for  each  node  and  link  within  the  superlink 

4)  Calculate  boundary  condition  coefficients 

5)  Using  all  of  the  above  coefficients  and  recurrence  relations,  calculate  a  set  of 
coefficients  for  use  in  the  solution  matrix 

6)  Repeat  steps  1-5  for  each  superlink  in  the  system 

7)  Calculate  matrix  values  based  on  connectivity  of  the  superlinks 

8)  Solve  the  sparse  matrix 

9)  Based  on  the  matrix  solution  of  head  at  each  unknown  junction,  calculate  the  flow 
at  the  upper  and  lower  ends  of  each  superlink 
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10)  Solve  for  flow  and  depth  in  each  pipe  and  node,  repectively 

11)  Continue  from  step  1  with  new  flows  and  depths 

Code  testing  and  refining 

Testing  of  the  scheme  began  with  the  simplest  network  possible;  one  superlink 
connecting  two  junctions  with  known  heads.  The  true  test  of  index  accuracy  was  to 
ensure  that  an  equivalent  flow  was  calculated  in  both  forward  and  reverse  directions  after 
simply  switching  the  boundary  heads.  Although  depth  is  determined  at  nodes,  area  of 
flow  is  required  at  the  center  of  the  links  to  calculate  flow.  A  depth  at  the  midpoint  of 
each  link  exists  solely  for  the  purpose  of  flow  calculations,  and  is  an  average  of  adjacent 
node  depths  (Figure  4).  Stability  problems,  particularly  those  related  to  flow  reversal, 
were  solved  by  revisiting  the  entrance/exit  hydraulics  as  described  in  the  following 
sections. 

Linki 

Node  1  Node  2 

Figure  4  SUPERLINK  depth  computation  scheme 

Entrance  hydraulics 

The  general  equation  for  inlet-controlled  flow  is  given  as 

Q  =  CA^2gNH  (6) 

where  C  is  a  geometric  coefficient,  A  is  the  flow  area,  and  NH  is  the  difference  in  head 
between  the  supply  reservoir  (junction)  and  pipe  (node  1,  link  1).  Ji  (1998)  had  taken 
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entrance  boundary  equations  from  other  sources,  and  thus  the  derivation  could  not  be 
easily  followed.  Re-deriving  the  superlink  end  equations  from  Eq.  5  created  an  alternate 
set  of  boundary  conditions.  We  define  A//  =  H  -h-  where  H  is  the  junction  head, 

h  is  the  depth  at  the  first  node,  and  Z^v  is  the  invert  elevation  of  the  first  node.  By 
squaring  both  sides  of  the  flow  equation  we  get; 

Q^=C^A-g(H-h-Zj  (7) 


The  time  varying  Q  is  broken  into  the  current  time  step  and  the  future  time  step,  and  we 
solve  for  depth  h.  (Eq.  8)  This  process  can  be  applied  to  the  downstream  end  of  a  pipe  as 
well  to  account  for  instances  of  backward  flow.  (Eq.  9) 
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With  these  changes  implemented  in  the  code,  the  model  ran  with  a  very  simple  one-pipe 
network  with  known  heads  at  either  end  junction. 


Exit  hydraulics 

Eike  the  pipe  entrances,  pipe  exits  were  modified  from  Ji’s  (1998)  algorithm  to 
more  accurately  model  various  flow  regimes.  In  exit  hydraulics,  four  possible  conditions 
must  be  considered  for  pipes  flowing  less  than  full.  Eirst,  in  a  mild  sloped  channel  where 
normal  depth  is  greater  than  critical  depth  and  junction  head  less  than  critical  depth,  the 
outfall  depth  is  controlled  by  critical  depth.  Second,  steep  sloped  channel  exits  where 
critical  depth  is  greater  than  normal  depth  are  governed  at  the  exit  by  normal  depth. 
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Third  is  a  critical  sloped  channel  where  normal  depth  equals  critical  depth,  critical  depth 
is  used  for  end  control.  Fourth  is  the  case  of  backwater  effects,  when  the  depth  in  the 
downstream  junction  begins  to  affect  depths  upstream.  To  experience  backwater  effects, 
the  head  in  the  junction  must  exceed  the  head  of  critical  depth  (critical  depth  +  invert 


elevation). 


Solid  line:  Mild  slope 
Dashed  line:  Steep  slope 


Figure  5  Pipe  exit  conditions 


If  the  system  is  obeying  conservation  of  momentum  and  the  length  of  pipe  is  sufficient 
such  that  the  friction  slope  is  equal  to  the  bed  slope,  the  solved  depth  should  be  normal 
depth.  Critical  depth,  however,  must  be  calculated  for  the  given  flow  rate  and  geometric 
variables.  As  the  solution  for  critical  depth  is  non-linear,  the  Newton-Raphson  iterative 
solution  technique  was  employed.  This  method  searches  for  roots  of  an  equation  with  a 
truncated  Taylor  series  expansion  to  approximate  F(x)  for  some  guess  x.  In  this  case,  the 
function  F(x)  is  Manning’s  equation  for  open  channel  flow.  (Eq.  10)  Because  the  series 
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is  truncated,  the  solution  is  not  perfect.  A  correction  is  applied  by  Equation  1 1 ,  where  Ax 


is  a  eorreetion  to  flowrate  Q. 


1  2/  1/ 

Q  =  -AR^^S'^\  (10) 
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/ 9x 


Iterations  are  performed  until  the  eorreetion  Ax  drops  below  a  given  threshold.  (Mays, 
1999).  This  proeess  was  imbedded  into  the  eode,  and  at  eaeh  time-step  the  exit  depth  is 
eheeked  against  both  eritieal  depth  and  downstream  junction  depth. 


Model  Verification 

To  verify  that  the  eomplete  model  was  operating  properly,  Ji’s  (1998)  test 
simulation  for  a  simple  six  pipe  network  was  reprodueed.  A  boundary  head  at  junetions 
D  and  F  represented  a  tidal  surge  indueing  baekwater  effeets.  Input  at  junetions  A  and  C 
began  after  the  tidal  surge,  allowing  negative  flow  to  oeeur.  Positive  pipe  slope  is 
denoted  by  the  small  arrows  along  eaeh  superlink. 


27 


The  following  hydrographs  matched  those  provided  by  Ji  (1998).  Reverse  flow  was 
present  in  all  links  at  the  start  of  the  simulation  due  to  the  downstream  boundary 
condition.  To  demonstrate  the  importance  of  looped  network  capability,  pipe / flowed 
opposite  its  slope  direction  throughout  the  simulation  because  of  a  greater  head  at 

junction  C  than  E. 
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Figure  7  Superlink  demonstration  simulation 


Conduita 
Conduit  b 


Changes  to  Model 

As  complexity  was  added  to  the  test  networks,  additional  deficiencies  in  the 


original  formulation  were  discovered.  One  such  instance  involved  the  introduction  of 
unknown  heads  to  the  model,  thus  incurring  the  use  of  the  matrix  solution.  A  mass 
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balance  check  indicated  model  connectivity  problems,  but  once  again  it  was  extremely 
difficult  to  track  the  source  of  error.  By  process  of  elimination,  it  was  discovered  that  the 
signs  of  two  summations  in  a  coefficient  equation  had  been  inadvertently  switched  in  Ji’s 
(1998)  publication.  In  similar  fashion  another  typographic  error  was  found  in  an  equation 
relating  boundary  conditions  where  a  parenthesis  had  been  omitted. 

The  last  major  hurdle  dealt  with  initial  conditions.  The  nature  of  equations  8  &  9 
does  not  allow  flow  to  move  into  the  system  when  the  area  of  flow  is  zero.  It  is  therefore 
necessary  to  maintain  a  very  small  depth  at  the  nodes  even  when  flow  is  zero.  A  danger 
in  imposing  a  depth  is  to  create  instability  within  the  flow  calculation,  as  physically  these 
numbers  should  be  simultaneously  generated.  Extensive  testing  found  that  an  initial 
depth  of  0.00001  m  provided  a  stable  minimum,  allowing  flow  to  commence  without 
disrupting  the  mass  and  energy  balance.  This  value  is  likewise  imposed  when  inputs 
cease  and  a  network  drains,  simply  to  keep  the  pipes  “wet”. 
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GSSHA  and  SUPERLINK 


Linking  Models 

With  a  robust  set  of  options  comes  complexity,  hence  the  GSSHA  code  greatly 
encourages  new  subroutines  to  run  side  by  side  instead  of  attempting  to  imbed  the 
processes.  This  was  beneficial  at  the  development  stage  (as  testing  could  be  limited  to 
SUPERLINK  exclusively)  as  well  as  the  linking  stage  (see  Appendix  A). 

Linking  the  two  models  presented  a  challenge  in  simulating  hydrologic  processes 
on  a  wide  range  of  scales.  Typical  GSSHA  grid  sizes  are  on  the  order  of  10-30  meters, 
which  is  usually  suitable  for  modeling  catchment  scale  processes  such  as  infiltration, 
evaporation,  overland  flow,  and  channel  flow.  However,  it  became  apparent  that  grate 
inlets  function  by  micro-topography,  as  curbs  and  crowned  road  profdes  direct  flow  to 
the  catch  basins.  This  scale  of  hydrology  cannot  be  represented  on  even  a  10-meter  grid, 
so  it  was  necessary  to  develop  a  conceptual  relationship  between  GSSHA  and 
SUPERLINK  to  realistically  introduce  flow  to  the  storm  drainage  network. 

Grate  Hydraulics 

Some  simplifying  assumptions  were  necessary  to  simulate  the  intake  process 
feasibly.  It  is  reasonable  to  assume  that  most  inlet  points  are  located  at  depressions  or 
curbs  of  a  crowned  roadway.  Although  this  certainly  is  not  true  in  all  cases,  the  vast 
majority  of  the  catch  basins  fit  this  description.  GSSHA  considers  planar  sloping  grid 
cells  with  a  depth  of  ponded  water  for  the  case  of  overland  flow.  At  each  intake  point, 
SUPERLINK  is  given  the  number  of  grates  represented  by  a  node  (A  =  1  to  4).  The 
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potential  flow  intake  available  to  SUPERLINK  {qin)  is  then  taken  as  a  pereentage  of  total 
ponded  volume  (Fee//)  based  on  the  number  of  grates  per  time  step  {dt)  by  the  equation 


NaV 

cell 

dt 


(12) 


where  a  =  l/Nmax-  This  simplilieation  follows  the  reasoning  that  a  eell  with  four  grates 
would  be  eapable  of  intereepting  all  of  the  ponded  water. 


Surcharged  manholes  &  grates 

It  is  elear  that  the  storm  drainage  network  and  the  land  surfaee  seheme  must 
internet  in  sueh  a  manner  as  to  allow  sureharged  manholes  and  grates  to  diseharge  baek  to 
the  surfaee.  As  the  SUPERLINK  model  will  run  in  parallel  with  GSSHA,  this  interplay 
is  easily  aeeounted  for.  Eaeh  time  SUPERLINK  is  ealled,  the  depth  in  the  eells 
eontaining  manholes  and  grates  will  be  passed  to  the  subroutine.  SUPERLINK 
determines  whether  there  is  suffieient  spaee  in  the  strueture  to  aeeept  all  flow  available 
from  the  grate  hydraulies  ealeulation.  If  so,  this  volume  will  be  subtraeted  from  the 
overland  flow  plain  and  added  to  the  subsurfaee  network.  Otherwise,  the  flow  will  be 
foreed  to  eontinue  on  in  the  overland  plane,  thus  simulating  a  system  operating  at  full 
eapaeity.  If  SUPERLINK  ealeulates  the  head  at  a  manhole  to  be  greater  than  the  ground 
surfaee  elevation  or  overland  flow  head,  the  exeess  volume  will  be  subtraeted  from  the 
drainage  network  and  added  to  the  surfaee  eell.  By  following  this  method  at  eaeh 
timestep,  SUPERLINK  and  GSSHA  will  interaet  realistieally. 
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IV.  METHODOLOGY 


Data  for  the  Modeling  Study 

It  was  necessary  to  model  a  watershed  with  a  significant  urban  presence  and 
subterranean  drainage  network  to  fully  test  the  routines.  A  low  gradient  topography 
would  provide  situations  of  inundation  and  pressurized  pipes.  But  perhaps  most  critical 
was  the  availability  of  a  quality  dataset  including:  rainfall  records,  stream  flow  records, 
digital  elevation  model  (DEM),  land-use  coverage,  stream  channel  and  storm  drainage 
network  data.  All  but  one  of  these  requirements  could  be  readily  met  by  Dead  Run,  a 
14.3  sq  km  watershed  in  Baltimore,  Maryland. 


Figure  8  Maryland  with  location  of  14.3  km  Dead  Run  watershed 
Bounding  UTM  Coordinates:  4355809  -  4349732  N,  347633  -  352363  E 


Extensive  field  work  has  been  ongoing  as  part  of  a  study  by  researchers  from 
Princeton  University.  Precipitation  data  were  derived  from  the  WSR-88D  radar  in 
Sterling,  Virginia  (75  km  from  watershed  center)  and  a  network  of  19  rain  gages. 
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Discharge  data  at  the  outlet  of  Dead  Run  is  provided  by  a  USGS  gaging  station.  The 
Prineeton  group  has  obtained  GIS  layers  of  land  use,  a  10  meter  DEM,  and  drainage 
maps  from  the  eounty  of  Baltimore.  Their  own  field  eampaigns  have  involved  eross 
seetion  measurements  of  stream  ehannels,  rating  eurve  measurements,  and  soil  infiltration 
testing.  The  basin  also  ineludes  a  number  of  small  to  medium  sized  detention  basins. 
Only  one  element  required  a  signifieant  amount  of  work:  the  drainage  network. 

Storm  drainage  network 

The  three  maps  that  eovered  this  watershed  eontained  a  plan  view  of  the  network 
with  symbols  for  eatch  basins,  manholes,  outlets,  and  pipe  sizing.  To  produee  pipe  data 
for  use  in  the  SUPERLINK  seheme,  it  was  neeessary  to  sean  and  georeferenee  these 
images.  Eollowing  is  an  example  of  the  seanned  image  and  resulting  digitized  network 
superimposed  on  the  GIS  land  use  map.  The  GIS  landuse  map  (Eigure  10)  also  shows  the 
remarkable  level  of  detail  available  for  this  watershed. 
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Figure  9  Scanned  drainage  map 


Figure  10  Digitized  drainage  map 


The  storm  drain  maps  were  created  independently  of  the  current  GIS  road  and 
building  layers,  and  did  not  overlay  each  other  perfectly.  This  is  due  to  the  fact  that  the 
pipe  network  scans  were  from  paper  maps  created  by  hand  and  not  tied  to  a  coordinate 
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system.  For  overland  flow  in  GSSHA  and  pipe  flow  in  SUPERLINK  to  internet 
realistieally,  network  intake  points  (eateh  basins)  should  coineide  with  impervious 
roadway  eoverage  where  applieable.  Thus  the  loeation  of  speoifie  network  elements  had 
to  be  shifted  slightly  to  aeeurately  represent  the  existing  eonditions.  Junetions  are  shown 
with  star  symbols,  nodes  (eateh  basins  and  eomputational  points)  by  open  eireles,  and 
superlinks  by  dashed  lines. 

The  numbering  seheme  required  by  SUPERLINK  is  diseussed  in  the  appendix. 
Junetions,  links,  and  nodes  were  ereated  as  GIS  layers  in  AreView  to  represent  pipe 
interseetions,  pipes,  eateh  basins,  and  outlets.  Beeause  this  was  the  first  use  of  the  storm 
drainage  model,  all  dataset  preparation  was  done  manually  to  determine  the  exaet 
ereation  proeess  of  the  input  file.  Eor  this  seenario,  all  ealeulations  were  performed  using 
points  (junetions  and  nodes)  with  assoeiated  tables  listing  downstream  pipe  size. 
Conneetivity  information  was  eontained  within  line  layers  that  simply  denoted  the 
upstream  and  downstream  junetion  number.  Einal  assimilation  of  this  GIS  data  is  further 
deseribed  in  Appendix  B. 

Pipe  invert  elevations  were  not  readily  aeeessible  for  the  Dead  Run  eatehment 
exeept  in  detailed  engineering  drawings.  Due  to  the  amount  of  data  being  proeessed,  a 
simplifying  assumption  was  used  to  ealeulate  the  invert  of  eaeh  pipe.  It  is  reasonable  to 
assume  that  a  typieal  subterranean  network  largely  follows  the  ground  surfaee  based  on 
pipe  installation  limitations.  With  this  in  mind,  AreView  was  used  to  assign  a  ground 
surfaee  elevation  to  eaeh  point  from  the  DEM.  Using  a  typieal  pipe  depth  of  2  meters, 
an  invert  was  ealeulated  from  the  surfaee  elevation.  Where  it  was  obvious  that  this 
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simplification  defied  engineering  praetices,  the  invert  elevations  were  modified  to 
produee  positive  slopes  in  the  assumed  flow  direetion. 

Hurricane  Isabel 

The  drainage  network  will  have  its  most  pronounced  effect  during  a  moderate  to 
high  intensity  storm,  and  therefore  the  event  of  September  18-19,  2003  was  seleeted  as 
the  model  test  case.  Hurricane  Isabel  was  elassified  as  a  Category  2  when  it  struck  the 
North  Carolina  coast,  but  weakened  to  a  tropieal  depression  as  it  moved  north.  The  storm 
dropped  heavy  rainfall  on  mueh  of  the  east  eoast,  including  Maryland  and  the  Baltimore 
watershed  of  Dead  Run.  As  is  eommon  in  hurricane  precipitation  patterns,  Baltimore 
received  two  strong  pulses  of  rainfall  150  minutes  apart  (Figure  1 1).  The  peak  discharge 
recorded  by  the  USGS  gaging  station  at  the  outlet  of  the  watershed  was  just  under  40  ems 
(-1400  efs  or  10  mm/hr),  and  significant  flooding  was  reported  on  the  North  Branch  of 
the  Potomae  River.  Basin  averaged  rainfall  peaked  at  53  mm/hr,  but  loealized  eells  of 
intense  precipitation  were  recorded  above  200  mm/hr.  Thus  the  distributed  nature  of  the 
rainfall  input  is  as  eritieal  as  the  distributed  landuse  and  soil  classifieation.  This  event 
was  seleeted  because  it  allows  ealibration  of  GSSHA  using  the  observed  preeipitation  and 
stream  reeords. 

A  brief  analysis  of  the  volume  of  rainfall  versus  volume  of  discharge  in  each  peak 
highlights  a  problem.  There  is  25%  less  precipitation  in  the  second  pulse  versus  the  first, 
but  only  9.5%  less  runoff  volume  for  the  seeond  peak  of  the  hydrograph.  If  the  initial 
conditions  were  unknown,  this  imbalanee  could  be  attributed  to  dry  soil  eonditions 
leading  up  the  first  peak.  However,  the  steady  period  of  light  rain  in  the  hours  preceding 
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the  storm  reduces  this  possibility.  The  model  results  will  demonstrate  this  apparent 
deficit  in  greater  detail  and  point  to  radar  bias  as  the  problem. 

60 


Figure  12  provides  a  sense  of  rainfall  intensity  compared  to  the  infiltration 
capacity  of  the  soil.  Rainfall  records  prior  to  the  simulation  period  are  similar  to  the 
period  from  0  to  150  minutes,  or  nearly  equal  to  the  saturated  hydraulic  conductivity. 

This  observation  is  useful  in  initializing  soil  moisture.  The  longer  first  pulse  of  rainfall 
produced  moderate  runoff,  while  the  shorter  second  pulse  generated  somewhat  more 
intense  infiltration  excess.  Although  the  storm  had  periods  of  remarkable  intensity,  storm 
total  accumulation  was  only  66  mm  (2.6  in)  and  therefore  is  not  categorized  as  an 
extreme  event. 
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Figure  12  Basin  average  rainfall  normalized  by  saturated  hydraulic  conductivity 
As  previously  mentioned,  radar  data  from  Sterling,  Virginia  was  provided  at  1  km 
spatial  and  6  minute  temporal  resolution,  and  overlays  the  Dead  Run  watershed  as  shown. 


Basin  outlet 


Figure  13  1  km  radar  grid  overlaying  watershed  (south-east  corner  of  each  cell  shown) 
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Parameter  Assignment 


Gridded  Watershed 

With  the  modeling  data  in  various  formats,  the  last  step  is  to  assemble  the 
components  for  use  within  the  modeling  framework  by  means  of  the  aforementioned 
Watershed  Modeling  System,  or  WMS.  This  software  is  capable  of  transforming 
polygon,  line,  and  point  data  into  gridded  maps,  which  are  then  read  by  the  GSSHA 
routine.  Initial  efforts  were  focused  on  a  10  meter  grid  size,  as  this  matches  the 
resolution  of  the  DEM  from  Princeton.  However,  each  hydrologic  process  must  be 
executed  on  the  resulting  143,000  grid  cells  (within  the  14.3  sq.  km  watershed),  and  this 
is  computationally  expensive.  Also,  gridded  models  are  subject  to  artifacts  in  the  DEM 
that  do  not  allow  the  land  surface  to  drain  properly.  Removal  of  this  “DEM  noise”  is 
much  more  difficult  at  10  meters  compared  to  30  meters,  because  the  larger  grid  size  is 
smoothed  during  aggregation.  A  30  m  grid  size  reduces  the  number  of  grid  cells  to 
approximately  16,000,  is  free  of  artifacts,  and  thus  was  chosen  for  the  initial  watershed 
simulations. 

A  landuse  classification  index  is  then  assigned  to  each  grid  cell  based  on  the 
majority  coverage  overlaying  it.  Tied  to  each  index  value  are  hydrologic  parameters  such 
as  saturated  hydraulic  conductivity,  roughness,  porosity,  initial  moisture  content,  and 
capillary  head. 

Table  1  Summary  of  Dead  Run  landuse 


Eand  Use  Classification 

2 

Area  (km  ) 

Percent  of 
Watershed 

Road  and  Parking  Eot 

1.57 

11.0% 

Building  Rooftop 

3.54 

24.8% 

Grass/Woodland 

9.19 

64.2% 
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Figure  14  30  meter  gridded  landuse  map 


Figure  15  GIS  landuse  map 
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The  gridded  landuse  index  map  (Figure  14)  retains  the  major  features  of  the 
watershed  seen  in  the  GIS  map  (Figure  15).  Even  at  the  30-meter  seale,  the  highway  that 
divides  the  watershed  into  quadrants  is  still  elearly  visible  as  well  as  the  large  buildings 
and  parking  lots.  Somewhat  less  distinetive  are  seeondary  roads  and  residential 
dwellings.  Despite  the  loss  of  resolution  due  to  eoarsening,  analysis  of  the  relative 
influenee  of  model  eomponents  will  still  be  aeeurate. 

Width  Function 

The  effeets  of  drainage  density  have  long  been  hypothesized,  but  the  powerful 
distributed  nature  of  this  model  provides  an  opportunity  to  simulate  various  degrees  of 
modifieations  and  eompare  their  relative  hydrologie  response.  To  quantitatively  report 
the  degree  to  whieh  a  watershed  is  drained  is  not  trivial  task,  beeause  the  measure  must 
take  into  aeeount  not  only  the  absolute  number  of  flow  links  but  also  their  length  to  the 
outlet.  For  example,  a  watershed  may  have  a  dense  network  of  ehannels,  but  if  arranged 
in  a  sinuous  fashion  will  respond  very  differently  to  a  basin  with  the  same  drainage 
density  but  straight  ehannels.  Thus  the  width  funetion  is  employed  to  develop  a  measure 
of  both  drainage  eharaeteristies.  For  this  researeh,  the  width  funetion  is  defined  as  the 
number  of  flow  segments  at  a  given  distanee  from  the  watershed  outlet.  A  flow  segment 
is  defined  as  a  10  m  length  of  either  ehannel  of  storm  sewer.  The  distanee  is  measured 
along  the  ehannels  and  sewers  themselves,  thus  representing  the  length  a  drop  of  water 
must  flow  to  the  exit. 

Three  seenarios  have  been  eonsidered  to  explore  the  effeet  of  the  width  funetion: 
natural  eonditions,  the  modified  ehannel  network,  and  the  eurrent  ehannel  network  with 
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storm  sewers.  To  simulate  the  pre-development  channel  configuration,  the  WMS  routine 
TOPAZ  (Garbrecht  and  Marz,  1993)  was  used  to  calculate  fiow  accumulations.  This 
algorithm  simply  determines  the  surface  area  that  drains  through  any  given  cell.  Given  a 
threshold  value,  WMS  then  creates  stream  arcs  from  cells  with  a  contributing  drainage 
area  above  this  minimum.  This  approximation  of  natural  conditions  is  shown  in  the 
following  Figure  16,  overlaid  on  the  complete  existing  network.  The  second  and  third 
cases  are  shown  in  Figure  17,  with  channels  in  blue  and  storm  sewers  in  red. 


Figure  16  Idealized  natural  network  Figure  17  Existing  channel  network  in  bold, 

storm  drainage  network  in  thin  lines 


The  width  function,  plotted  against  the  distance  from  the  outlet  in  Figure  18,  provides  a 
tool  to  qualitatively  assess  the  modifications  of  the  basin  drainage.  Comparing  the 
natural  and  existing  channel  histograms,  little  deviation  occurs  until  approximately  4000 
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meters.  At  this  point,  the  modified  channels  display  a  significant  increase  in  the  drainage 
density.  A  visual  inspection  of  the  change  in  blue  lines  from  Figure  16  to  Figure  17 
explains  the  difference;  numerous  channels  were  installed  in  the  central  core  of  the 
watershed.  The  existing  channel  plot  also  extends  1000  meters  beyond  natural 
conditions,  due  to  an  apparent  extension  of  the  channels  farthest  from  the  basin  outlet. 
Figure  18  shows  that  expansion  of  the  open  channel  network  has  a  significant  impact  on 
the  width  function. 

30 


Existing  Channels  &  Drainage 
Existing  Channels 


Distance  from  Outlet  (m) 

Figure  18  Effect  of  drainage  network  on  width  function 

The  effect  of  storm  drainage  on  the  width  function  is  extremely  pronounced. 
Once  again,  notable  increases  over  the  existing  channels  can  be  matched  to  heavily 
sewered  areas  in  Figure  17.  The  strongest  impact  again  falls  beyond  4000  m,  and  not 
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coincidentally.  In  regions  where  the  ehannels  have  been  modified,  it  is  likely  evidenee  of 
intensive  land  development,  whieh  is  often  aeeompanied  by  the  installation  of 
subterranean  drainage.  By  inspeetion  of  the  drainage  map,  regions  in  the  south  and  west 
of  the  watershed  have  the  highest  density  of  sewers.  This  is  refieeted  in  the  width 
funetion,  with  a  drastic  increase  in  the  number  of  flow  segments  past  the  eorresponding 
length  of  5000  m.  However,  eonelusions  on  the  effeet  of  the  inerease  in  density  must 
inelude  mention  of  eonveyanee.  The  eonveyanee  K,  or  quantity  of  flow  a  given  segment 
ean  pass,  is  much  less  in  the  storm  drainage  pipes  than  in  natural  ehannels. 

1  2/ 

K  =  -AR/^  (12) 

n 

In  this  equation,  n  is  Manning’s  roughness  eoeffieient,  A  is  the  bank-full  flow  area,  and  R 
is  assumed  to  be  the  bank-full  flow  hydraulie  radius.  For  pipes,  eonveyanee  is  ealeulated 
at  full  flow  eonditions  as  well.  The  influenee  of  eonveyanee  ean  be  signifieant.  For 
example,  a  trapezoidal  ehannel  1  m  wide  has  a  eonveyanee  of  400,  eompared  to  a  0.45  m 
pipe’s  eonveyanee  of  1.5. 

By  weighting  the  width  funetion  by  eonveyanee,  the  segments  no  longer  have 
equal  eount  values.  Though  the  large  order  streams  near  the  outlet  may  be  few  in 
number,  they  will  have  a  mueh  larger  eonveyanee  sum  than  the  numerous  first  order 
streams.  Although  the  inerease  in  density  seen  in  Figure  17  represents  lower  eapaeity 
segments,  they  are  still  apparent  in  the  eonveyanee  weighted  plot. 
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Figure  19  Comparison  between  normal  and  conveyance  weighted  width  functions 


The  increase  in  drainage  density  may  also  be  offset  by  another  feature  that  appears  in  the 
width  function  plot.  The  majority  of  the  density  increase  occurs  beyond  a  distance  where 
the  natural  system  registered  flow  segments  (6000  m).  This  can  be  explained  by  the 
indirect  routing  of  small  elements  of  the  drainage  system.  Hence,  a  longer  distance  of 
travel  to  the  outlet  will  partially  suppress  the  response  time  effect  of  dense  drainage. 


Hydrologic  processes 

The  importance  of  various  hydrologic  processes  to  a  modeler  depends  on  the 


focus  of  the  modeling  effort.  Long-term  simulations  are  heavily  influenced  by 
groundwater  contributions  to  streamflow  and  climatic  factors  such  as  evapotranspiration 
losses.  Single  event  simulations  targeting  flood  peaks  and  timing  rely  mainly  on  the 
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short-term  processes  of  infiltration,  overland  flow  routing,  and  channel  routing. 
Additionally,  the  effect  of  lakes  and  hydraulic  structures  are  important  for  both  long  and 
short-term  simulations.  As  the  focus  of  this  study  is  the  effect  of  drainage  network 
modifications,  the  model  components  will  include  the  processes  fundamental  to  single 
events. 

Setting  hydrologic  parameters  in  GSSHA  utilizes  both  values  from  the  literature 
and  past  experience  with  the  model.  Saturated  hydraulic  conductivity,  roughness 
coefficients  of  landuse,  and  initial  soil  moisture  are  the  key  parameters  for  the  processes 
involved  in  a  short-term  simulation.  One  soil  type  was  used  throughout  the  watershed 
because  soils  data  were  not  readily  accessible.  Of  the  two  soil  parameters,  the  model  is 
most  sensitive  to  saturated  hydraulic  conductivity,  and  thus  it  becomes  the  calibration 
parameter.  As  for  the  other  soil  parameter,  because  the  simulations  were  begun  after  a 
period  of  steady  light  rain,  the  soil  surface  was  assumed  to  be  saturated. 

Roughness  coefficients  influence  the  travel  time  across  overland  flow  planes  and 
in  the  channels.  Their  values  were  estimated  from  past  experience  with  GSSHA 
modeling  as  well  as  published  values.  Overland  flow  Manning’s  roughness  coefficients 
were  set  at  0.400  for  the  overland  plain  and  0.050  for  roadways  and  parking  lots.  The 
disparity  between  these  two  values  illustrates  the  vast  difference  between  paved  and 
natural  surfaces.  In  contrast  to  the  roadways  and  parking  lots,  a  rooftop  area  roughness 
coefficient  of  0.800  simulates  the  delay  due  to  gutters  and  roof  drains.  Both  of  these 
surfaces  do  not  experience  infiltration,  but  the  runoff  response  from  roofs  is  substantially 
slower  than  pavement.  Once  runoff  reaches  the  channels,  a  typical  value  of  0.04  is  used 
for  channel  routing  based  on  a  clean,  winding  stream  (Chow,  1988). 
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Model  Calibration 


To  validate  the  joint  model  formulation,  the  GSSHA  outflow  hydrograph  was 
compared  to  the  USGS  observed  discharge  record  at  the  basin  outlet.  Saturated  hydraulic 
conductivity  was  varied  between  0.1  and  0.8  cm/hr  as  simulation  results  were  compared 
to  the  observed  plot.  Figure  20  displays  the  final  simulation  with  a  Ksat  value  of  0.5 
cm/hr.  The  model  is  most  sensitive  to  Kgat,  and  the  final  value  was  chosen  for  its  ability 
to  produce  a  first  peak  of  40  cms.  Increasing  Ksat  to  0.8  cm/hr  reduce  this  first  peak  by 
almost  25%. 
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Figure  20  Calibration  of  GSSHA  on  Dead  Run 


It  is  evident  that  this  “manual  calibration”  did  not  produce  a  perfect  match  to  the 
observed  hydrograph.  A  system  such  as  Shuffled  Complex  Evolution  would  give  a  better 
parameter  set  by  simulating  thousands  of  combinations  of  variables  within  a  set  range. 
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But  parameters  are  not  the  only  reason  for  the  inability  of  GSSHA  to  mateh  the 
hydrograph,  and  more  speeifieally  the  seeond  peak.  Reealling  the  diseussing  of 
preeipitation  volume  versus  runoff  volume  from  the  model  data  seetion,  the  volume  of 
rainfall  in  the  seeond  pulse  was  signifieantly  less  that  the  first  pulse  while  their  runoff 
volumes  were  almost  the  same.  The  most  likely  explanation  of  this  is  radar 
underestimation.  Radar  data  are  almost  always  eorrelated  to  gage  data  by  a  bias 
adjustment.  A  mean  field  bias  between  radar  and  gage  data  was  used  by  the  Prineeton 
group  to  adjust  the  radar-derived  intensities  by  the  following  equation  =1.5  \r^  where 

Ta  is  the  adjusted  rate  and  ro  is  the  observed  rate.  Also,  radar  refieetivity  level  is  eapped 
at  speeified  deeibel  value  (55  dBZ)  to  prevent  highly  refleetive  hail  from  giving 
unreasonably  intense  rainfall  rates.  In  the  ease  of  Hurrieane  Isabel,  it  is  highly  probable 
that  a  eombination  of  these  two  faetors  produeed  a  radar  dataset  without  enough  rainfall 
during  the  seeond  phase  of  the  storm.  Without  question  the  highest  intensities  oeeurred 
during  the  latter  pulse,  thus  subjeeting  these  rates  to  the  potential  of  being  elipped  and  or 
dampened  by  the  mean  field  bias.  Ogden  et  al.  (2000)  showed  that  applying  one  radar 
bias  value  to  a  multi-pulse  storm  does  not  agree  well  with  gage  data  for  individual  pulses. 

The  four  detention  ponds  simulated  were  based  on  rough  approximations  of  size, 
outlet  strueture,  and  diseharge  eurves.  A  better  understanding  of  the  ponds  and  their 
attenuation  oharaeteristies  might  signifieantly  improve  modeling  results.  In  similar 
fashion,  a  thorough  geometry  of  the  ehannels  was  not  known,  thus  introdueing  another 
possibility  for  model  misrepresentation.  Shorteomings  due  to  data  in  the  model 
ealibration  are  not  as  eritieal  as  if  the  watershed  of  interest  were  being  simulated  for 
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design  purposes.  Rather,  its  utility  lies  in  the  ability  to  simulate  one  eatehment  with 
varying  drainage  densities  and  properties. 
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V.  RESULTS  &  DISCUSSION 


Storm  Drainage  Network 

With  the  model  development  eomplete,  it  is  now  possible  to  step  through  an 
idealized  timeline  of  the  urbanization  of  the  Dead  Run  watershed.  From  the  estimated 
natural  eonditions  to  eurrent  developed  and  drained  eonditions,  it  is  possible  to  simulate 
the  effeet  of  separate  and  joint  modifieations  to  the  basin.  Outlet  hydrographs  will 
demonstrate  the  effeets  of  man-made  hydrologie  features  ineluding  ehannel 
modifieations,  detention  basins  and  eulverts,  storm  drains,  and  impervious  area. 

Channel  Modifications  and  Detention  Basins 

Simulating  the  pre-development  watershed  is  elearly  a  subjeetive  task.  It  is 
impossible  to  reproduee  every  hydrologie  feature  that  existed,  say,  200  years  ago.  There 
are  some  approximations,  however,  that  ean  be  made  based  on  typieal  development 
trends  and  existing  topography.  As  mentioned  earlier,  the  ehannel  network  for  this 
simulation  was  derived  from  a  flow  aeeumulation  algorithm.  The  width  and  depth  of 
these  “natural”  ehannels  were  assigned  identieal  values  to  the  modified  network  to  avoid 
storage-related  differenees.  The  one  property  given  a  slightly  different  value  was  the 
roughness  eoeffieient.  Urbanized  ehannel  networks  often  are  intensively  maintained, 
redueing  flow  attenuating  material  sueh  as  brush  and  long  grass  whieh  might  persist  in  a 
natural  ehannel.  They  also  tend  to  be  straighter  and  more  effieient  than  a  naturally 
formed  ehannel.  Thus  a  Manning’s  n  of  0.05  was  used  rather  than  0.04,  as  reeommended 
by  Chow  (1988). 
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Three  networks  were  simulated  without  the  existing  impervious  areas  so  as  to 
remove  land  surface  factors  from  influencing  a  strictly  channel  comparison.  Figure  21 
shows  three  conditions:  “natural”  channels,  expanded  channel  network  as  exists 
currently,  and  current  channel  network  with  existing  detention  basins  and  culverts. 

40  ■  -  — 

- Detention  Basins  and  Culverts 

- Modified  Channels 

-  Natural  Conditions 


Figure  21  Effect  of  modified  channels  on  non-imp  envious  watershed 


Most  striking  about  Figure  21  is  the  change  from  a  pre-development  watershed  to 
modified  channels.  The  first  flood  peak  increases  by  67%,  as  runoff  volume  decrease 
minimally  (2.5%).  Comparing  the  thick  line  channels  in  Figures  16  and  17,  as  well  as 
their  corresponding  width  functions  in  Figure  18  provides  clear  evidence  that  increases  in 
main  channel  density  has  a  profound  effect  on  the  hydrograph.  The  addition  of  lakes  and 
culverts  reduces  the  peak  discharge  to  a  50%  increase  over  the  natural  case.  This  behaves 
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as  expected  due  to  the  storage  capacity  of  these  small  detention  ponds  and  the  attenuating 
effect  of  road  culverts. 


Table  2.  Model  results  of  channel  modifications,  no  imperviousness 


Scenario 

Peak  Discharge 
(cms) 

Time  to  First 
Peak  (min) 

Discharge 
Volume  (m^3) 

Natural  Channels 

19.3 

295 

322,000 

Modified  Channels 

32.3 

266 

314,000 

Lakes  &  Culverts 

28.9 

270 

303,000 

Time  to  the  first  peak  is  decreased  by  the  modifications  to  the  channel  network  (Table  2). 
Because  there  is  no  impervious  coverage  in  these  scenarios,  runoff  from  regions  of  the 
watershed  without  any  channels  is  significantly  delayed.  Despite  the  fact  that  no 
imperviousness  exists  in  these  three  scenarios,  increased  drainage  density  is  able  to 
convey  runoff  quicker  than  overland  flow. 


Storm  Drainage  Network 

The  next  phase  of  simulations  includes  the  existing  network  of  modified  channels, 
distributed  impervious  coverage,  subsurface  drainage,  and  detention  ponds.  Of  particular 
interest  is  the  performance  of  the  storm  sewer  component,  and  ultimately  its  effect  on  the 
hydrology  of  Dead  Run.  A  subset  of  the  268  links  has  been  selected  to  display  typical 
pipe  flow  hydrographs  (Figure  22).  The  majority  of  pipes  are  first  order,  receiving  all 
input  from  surface  grate  structures  (solid  lines).  The  increasing  length  of  dashes  on  the 
pipe  hydrograph  lines  represent  increasing  order  pipes,  accepting  inflow  from  both 
upstream  pipes  as  well  as  inlet  grates.  The  pipe  hydrographs  have  a  remarkably  quick 
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response  to  the  rainfall,  and  are  largely  drained  before  the  peak  discharge  reaches  the 
outlet. 
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Figure  22  Outlet  and  sample  pipe  discharge  hydrographs,  longer  dashes  represent 

higher  order  pipes 

The  basin  average  rainfall  hyetograph  has  been  provided  at  the  top  of  Figure  22 
for  comparison  with  pipe  flow.  The  principal  observation  from  this  is  the  pipe  flow 
during  the  period  of  intense  rainfall  during  the  second  pulse  at  ~350  min.  This  can  be 
explained  by  subsurface  hydraulics.  Once  pipes  become  full  and  flow  under  pressure, 
they  are  unable  to  accept  additional  input.  Also,  pressurized  manholes  and  intake  grates 
behave  as  outlet  points  when  the  head  exceeds  the  ground  surface  elevation.  Thus  this 
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Outlet  Discharge  (cms) 


period  of  extraordinarily  heavy  precipitation  (local  cells  up  to  200  mm/hr)  overwhelms 
the  network  and  much  of  the  runoff  remains  on  the  surface.  Although  the  first  storm 
pulse  also  produced  high  rain  rates,  the  storm  sewers  were  able  to  accept  most  of  the 
water  available  at  the  intake  grates. 

To  fully  explore  the  hypothesis  that  storm  sewers  become  overwhelmed  and  have 
less  of  an  effect  on  extreme  events,  the  1997  event  that  caused  flash  flooding  in  Fort 
Collins,  Colorado  (Ogden  et  ah,  2000)  was  simulated  over  the  Dead  Run  watershed.  This 
four-pulse  storm  dropped  14.98  cm  (5.9  in)  of  rain  in  4  V2  hours,  more  than  doubling  the 
flood  peak  of  the  Dead  Run  outlet  hydrograph  compared  to  Flurricane  Isabel.  Although 
the  draining  of  the  basin  was  slightly  enhanced  as  seen  in  Figure  23,  the  flood  magnitude 
was  unaffected.  The  limited  conveyance  and  intake  potential  of  the  drainage  network 
could  not  enhance  flow  routing  to  the  outlet  enough  to  raise  the  peak  discharge. 
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Figure  23  Reduced  effect  of  storm  sewers  with  extreme  event 
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Imperviousness 


With  the  drainage  network  fully  functional,  it  is  possible  to  explore  the  relative 
effect  of  various  model  components.  To  assess  the  question  of  impervious  coverage 
versus  storm  drainage,  Figure  24  contains  plots  of  sequentially  added  impervious  areas 
and  storm  sewers.  Case  1:  no  impervious  areas,  no  storm  sewer;  Case  2\  distributed 
impervious  areas,  no  storm  sewer;  Case  3:  distributed  impervious  areas  with  storm  sewer 
network.  All  cases  include  the  same  channel  network,  detention  basins,  and  culverts. 


Flood  magnitudes  are  amplified  by  the  addition  of  both  impervious  coverage  and 
the  storm  drainage  network.  The  first  peak  increased  9.1%  from  Case  1  to  Case  2,  and 
44%  from  Case  1  to  Case  3.  This  demonstrates  the  ability  of  the  storm  sewers  to 
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drastically  alter  the  peak  flood  hydrograph.  The  seeond  peak  does  not  react  as  much  to 
the  addition  of  storm  drains.  This  ean  be  attributed  to  the  previous  discussion  of  the 
seeond  pulse’s  preeipitation  intensity,  whieh  inundated  the  system.  Differenees  in 
increase  of  discharge  volume  are  not  as  pronounced.  From  Case  1  to  Cases  2  and  3,  the 
inereases  were  31%  and  39%  respeetively.  From  this  observation,  the  storm  drainage 
network  does  not  inerease  diseharge  volume  significantly  compared  to  impervious 
eoverage. 


Table  3  Results  from  modified  drainage  network  simulations 


Scenario 

Peak  Discharge 
(cms) 

Time  to  Peak 
(min) 

Discharge 
Volume  (m^3) 

No  imperviousness 

28.7 

270 

298,000 

With  imperviousness 

31.3 

268 

391,000 

Storm  Drainage 

41.2 

267 

416,000 

It  is  difficult  to  isolate  the  effects  of  impervious  coverage  and  storm  drainage 
networks.  They  are  not  independent  proeesses,  as  impervious  areas  feed  subsurface 
culverts  and  pipes.  Although  physically  unrealistic,  it  is  possible  to  model  them 
completely  separately.  Figure  25  shows  these  two  seenarios:  Case  T.  An  unsewered 
watershed  with  distributed  impervious  coverage  (isolates  storm  drainage);  Case  2\  A 
completely  pervious  watershed  with  the  storm  sewer  network  (isolates  imperviousness). 
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Figure  25  Isolated  effect  of  storm  drainage  compared  to  impervious  areas 


The  only  truly  useful  observation  from  this  set  of  simulations  oecurs  at  the  first 
peak.  Despite  the  lack  of  any  runoff  generating  imperviousness  in  Case  2,  the  drainage 
network  still  produces  a  higher  flood  magnitude  than  Case  1 .  The  fact  that 
imperviousness  and  the  corresponding  increased  drainage  density  both  play  an  integral 
role  is  not  disputed.  But  by  this  comparison,  it  is  evident  that  the  drainage  network  has 
the  capacity  to  produce  a  greater  effect  on  flood  peaks  than  imperviousness. 
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Drainage  Density 


TOPAZ  Networks 

Study  of  the  relationship  between  drainage  density  and  width  funetion  required 
developing  a  new  set  of  channel  networks.  While  using  the  existing  drainage  network 
would  have  been  preferable,  not  enough  stages  of  development  were  available  for 
meaningful  analysis.  Thus,  TOPAZ  was  once  again  employed  to  generate  6  increasingly 
dense  networks,  as  shown  in  Figure  26  (a)-(f).  To  explore  the  full  range  of  densities, 
channels  were  generated  by  3.0  km^  to  0.02  km^  accumulation  thresholds.  Comparing 
Figure  26(a)  to  Figure  16  shows  that  this  set  of  networks  begins  even  less  dense  than  the 
“natural”  network  used  for  previous  simulations. 


Figure  26  Increasingly  dense  drainage  networks 
Dd  =length  of  channel  /  watershed  area 
Aa  =accumulation  threshold 
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Figure  26  (cont) 
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Figure  27  Equally  weighted  width  functions  for  6  networks  without  storm  sewers 

The  most  valuable  observation  from  the  width  function  plot  concerns  the 
proximity  of  channels  to  the  outlet.  Drainage  density  increases  by  adding  additional 
channels  to  the  network.  Given  the  shape  of  the  Dead  Run  watershed,  these  new 
channels  occur  increasingly  farther  from  the  outlet.  In  other  words,  the  central  core  of 
the  watershed  is  already  almost  completely  drained,  and  few  links  are  added  in  this 
portion.  Thus  the  distribution  of  flow  segments  becomes  more  and  more  negatively 
skewed.  The  critical  statistical  measure  of  the  width  function  is  the  distance  to  the  mean. 
This  will  provide  a  scalar  term  describing  how  far  from  the  outlet  most  flow  segments 
occur.  From  visual  inspection,  a  less  dense  network  is  expected  to  have  a  smaller  mean 
than  a  very  dense  network. 
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Drainage  Density  Simulations 


Determining  the  effect  of  both  network  density  and  mean  width  function  requires 


a  GSSHA  simulation  for  each  test  case.  These  were  conducted  with  identical  distributed 


landuse  (with  impervious  areas),  roughness  values,  and  no  storm  sewers.  The  channel 


sizes  were  kept  equal  from  one  network  to  the  next  to  eliminate  the  influence  of 


conveyance  differences. 


Figure  28  GSSHA  simulation  hydrographs  for  each  drainage  density,  no  storm  sewers 
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Table  4  Results  from  density  simulations 


Accumulation 
Area  (km  ) 

Drainage 

Density  Dd 

2 

(km/km  ) 

Distance  to 
W(x)  Mean 
(m) 

Peak 

Discharge 

(cms) 

Time  to 
Peak  (min) 

Discharge 
Volume  (fn) 

0.02 

4.9 

4464 

65.0 

257 

445,000 

0.05 

2.8 

4048 

61.6 

257 

443,000 

0.10 

2.0 

3716 

57.6 

257 

437,000 

0.20 

1.5 

3353 

47.1 

256 

418,000 

0.50 

0.94 

2867 

32.6 

259 

399,000 

3.00 

0.39 

1714 

20.2 

259 

285,000 

Drainage  density  clearly  has  a  strong  effect  on  peak  discharge,  as  well  as  a  minor  impact 
on  overall  volume.  The  following  figure  displays  both  peak  discharge  and  volume  from 
each  of  the  six  density  scenarios.  The  drainage  density  term  has  units  of  kilometers  of 
channel  per  sq  km  of  watershed. 


Figure  29  Ejfect  of  channel  density  on  flood  peaks 
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2 

There  is  a  period  of  rapid  inerease  from  0.5  to  2  km/km  ,  but  then  levels  off 
beyond  3  km/km^  where  overland  flow  lengths  eonverge.  For  this  partieular  watershed, 
the  range  of  drainage  density  values  to  whieh  flood  peaks  are  very  sensitive  lies  less  than 
2  km/km  .  Beyond  this  range,  additional  density  does  not  influenee  flood  peaks  as 
strongly.  It  appears  that  that  flood  peaks  are  less  affeeted  by  drainage  network  expansion 
onee  developed  past  a  eritieal  density.  This  eould  be  partieularly  signifieant  for  a 
suburban  watershed  without  major  modifieations  to  the  natural  network.  If  its  drainage 
density  were  still  on  the  lower  end  of  the  sensitive  range,  relatively  minor  development 
eould  signifieantly  inerease  flood  magnitudes. 

Volumetrie  effeets  from  ehannel  density  are  not  as  pronouneed  (Figure  29). 
Although  densities  below  1  km/km^  display  a  eonsiderable  deerease  in  diseharge, 
inspeetion  of  the  hydrograph  (Figure  28)  shows  flow  rates  well  above  the  other  densities. 
Sinee  a  600  minute  simulation  time  did  not  allow  for  complete  draining  of  the  watershed 
at  these  very  low  densities,  conclusions  should  not  be  drawn  from  them.  For  the  more 
reasonable  densities  above  1  km/km  ,  runoff  volume  increases  only  slightly.  As  the 
drainage  network  expands  and  becomes  more  efficient  at  intercepting  overland  flow  and 
transporting  it  to  the  outlet,  there  is  reduced  infiltration  opportunity.  This  volume 
becomes  part  of  the  heightened  flood  peak  as  seen  in  Figure  28. 
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Width  Function 


Uniform  Density  Distribution 

As  previously  defined,  distanee  to  the  width  funetion  mean  is  a  deseription  of 
spatial  distribution  of  links.  The  TOPAZ  generated  networks  used  throughout  this 
analysis  were  generated  based  on  eontributing  area,  and  this  inherently  produees  uniform 
density  aeross  the  entire  watershed. 


Distance  to  Width  Function  Mean  (m) 


Figure  30  Effect  of proximity  of flow  segments  to  outlet  on  flood  peaks 


At  first  glanee.  Figure  30  appears  to  prove  that  flood  peaks  inerease  as  the  mass  of  flow 
ares  move  away  from  the  eatehment  outlet.  However,  inereasing  density  uniformly  over 
the  basin  forees  the  width  funetion  mean  to  move  further  out.  Therefore,  the  inerease  in 
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flood  peaks  is  due  to  the  drainage  density,  and  not  the  mean  distanee.  Conclusions  about 
the  effect  of  spatial  distribution  of  links  cannot  be  drawn  from  these  simulations. 

Distributed  Drainage  Density 

Two  scenarios  with  identical  drainage  density  must  be  considered  to  fully  explore 
the  effect  of  non-uniform  development  within  a  catchment.  Beginning  with  the  densest 
case  (Figure  26  f),  channels  were  removed  from  either  the  outer  or  central  regions  of  the 
watershed.  Selecting  certain  regions  to  remain  dense  produced  these  two  cases  (Figure 
31).  With  the  drainage  density  Dd  equal  in  both  cases,  the  effect  of  spatial  variability  is 
simulated  without  a  density  bias. 
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Figure  32  Width  function  plots  from  spatial  extremes 


The  width  functions  displayed  in  Figure  32  appear  as  expected,  as  each  scenario 
dominates  a  portion  of  the  plot  based  on  proximity  to  the  outlet.  It  is  clear  from  this  that 
identical  densities  can  produce  very  different  mean  distances.  The  normal  distribution 
for  density  close  to  the  outlet  has  a  mean  of  2864  m,  but  the  highly  skewed  distribution 
for  density  far  from  the  outlet  yields  a  mean  of  4043  m.  Given  the  drastic  differences  in 
mean  width  function,  a  significant  impact  on  the  hydrology  and  outlet  hydrograph  would 
be  expected. 

The  simulations  displayed  in  Figure  33  contradict  this  hypothesis.  In  fact,  the 
first  peak  is  virtually  unaffected  by  the  radical  change  in  the  drainage  network.  Although 
the  second  peak  does  exhibit  a  faster  rise  and  slightly  higher  peak,  the  effect  is  much  less 
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than  anticipated.  Less  travel  time  in  the  channel  for  the  ease  of  density  close  to  the  outlet 
allows  runoff  from  the  intense  second  pulse  of  rainfall  to  reach  the  outlet  quicker,  but 
substantial  volume  effects  are  not  evident.  The  eonclusions  from  these  results  are  simple: 
distribution  of  development  within  this  particular  watershed  does  not  seem  to  have  a 
pronouneed  effeet  of  flood  magnitudes,  and  eloser  drainage  densities  to  the  outlet  ean 
slightly  reduce  the  time  to  peak.  It  must  be  noted  that  the  realitely  small  size  of  the  Dead 
Run  watershed  may  be  eritieal  to  this  conclusion.  Large  watersheds  in  whieh  channel 
travel  time  plays  a  larger  role  may  show  very  different  response  to  spatial  distribution  of 


drainage. 


Figure  33  Ejfect  of  density  spatial  distribution  on  watershed  with  impervious  areas 
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Running  identical  simulations  on  a  watershed  without  any  impervious  coverage 
tested  the  influence  of  impervious  area.  Thus  four  cases  exist:  Case  1 :  Density  far  from 
outlet,  with  impervious  areas,  Case  2:  Density  near  the  outlet,  with  impervious  areas. 
Case  3:  Density  far  from  outlet,  without  impervious  areas.  Case  4:  Density  near  outlet, 
without  impervious  areas. 


Figure  34  Ejfect  of  density  spatial  distribution  on  watershed  without  impervious  areas 

Comparing  Figures  33  and  34  demonstrates  that  distributed  impervious  area 
reduces  the  effect  of  spatially  varied  drainage  density.  For  the  Cases  3  &  4,  without 
impervious  area,  significant  differences  are  evident  between  the  two  density  variations, 
whereas  little  contrast  was  apparent  in  the  Cases  1  &  2.  Since  much  of  the  impervious 
area  is  located  at  the  extremes  of  the  catchment,  streets  and  parking  lots  in  Case  2  could 
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be  reducing  overland  flow  times,  enabling  a  close  match  to  Case  1 .  For  the  Cases  2  &  4, 
with  density  near  the  outlet,  the  decreases  in  lag  time  are  approximately  equal.  But 
peaks,  especially  the  first,  are  sharply  increased  by  moving  the  drainage  density  farther 
from  the  outlet  in  Case  3 .  Shorter  average  overland  flow  lengths  will  allow  less 
infiltration,  thus  increasing  the  total  volume  of  discharge.  The  increased  channel  travel 
length,  however,  slightly  delays  this  increased  flood  peak. 

Storm  Drainage  to  Accumulation  Comparison 

Because  of  the  lack  of  storm  sewer  modules  in  many  current  distributed 
physically  based  models,  subterranean  drainage  pipes  are  often  approximated  by  open 
channels.  By  comparing  the  densest  TOPAZ  network  (Figure  35  b)  to  the  existing  Dead 
Run  network  with  storm  drainage  (Figure  35  a),  the  validity  of  this  approximation  can  be 
tested. 
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(a)  Dd  =  5.5  km/km^  (b)  Dd  =  4.9  km/km^ 

Figure  35  Drainage  network  to  flow  accumulation  comparison 

The  drainage  density  Dd  =  5.5  km/km  for  the  existing  storm  drainage  network  is 
appreeiably  more  dense  than  the  0.02  km^  aeeumulation  threshold  network.  Figure  36 
exhibits  the  similarities  of  the  two  width  functions.  The  distance  to  the  mean  width 
function  is  likewise  greater  for  the  existing  system,  at  5098  m  versus  4464  m.  By  the 
preceding  arguments,  it  would  follow  that  the  existing  system  should  produce  higher 
peaks. 
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Figure  36  Width  function  of  drainage  network  compared  to  flow  accumulations 
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However,  the  channelized  network  resulted  in  flood  peaks  over  65  cms,  whereas 
the  existing  network  simulation  produced  merely  41  cms  at  its  peak.  This  disparity  in 
peak  flow  values  is  evidence  of  the  profound  differences  between  natural  open  channels 
and  subterranean  pipes.  The  explanation  of  this  dissimilarity  comes  in  two  parts.  First, 
lateral  inflow  is  accepted  along  the  entire  length  of  an  open  channel  but  limited  to  inlet 
grates  for  a  storm  sewer.  Second,  conveyance  in  a  pipe  is  far  less  than  even  a  small  open 
channel  because  a  pipe’s  enclosed  geometry  limits  high  flows.  It  is  clear  through  the 
flood  magnitude  discrepancies  that  modeling  storm  sewers  with  open  channels  is 
generally  not  a  sound  approximation. 
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Initial  Soil  Moisture 


Testing  the  hypothesis  that  antecedent  soil  moisture  has  a  significant  impact  on 
flood  response  was  performed  on  both  distributed  impervious  land  use  and  non 
impervious  land  use.  Simulations  results  with  initial  soil  degree  of  saturations  of  20%, 
60%,  and  100%  are  shown  in  Figures  37  &  38. 


Figure  37  Effect  of  antecedent  moisture  on  watershed  without  impervious  area 


The  effect  of  initial  soil  moisture  is  clear.  It  is  expected  that  the  effect  would  be 
greater  in  a  watershed  without  impervious  area  because  of  the  greater  influence  of  soil 
properties  (such  as  infiltration)  in  such  a  basin.  But  even  in  the  simulation  with 
impervious  areas,  peaks  increased  over  40%.  Inspection  of  the  peaks  in  Figures  37  and 
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38  show  that  the  difference  between  the  peaks  of  the  two  scenarios  lessens  as  the  initial 
conditions  approach  saturation  (100%).  A  fully  saturated  watershed  without  impervious 
area  begins  to  behave  much  like  one  with  impervious  area,  but  the  volume  removed  by 
infiltration  still  decreases  peak  flows. 


Figure  38  Ejfect  of  antecedent  moisture  on  distributed  watershed  with  impervious  areas 

The  effect  of  antecedent  moisture  is  very  pronounced  for  this  storm  because 
infiltration  still  plays  a  significant  role.  Flowever,  the  influence  of  initial  soil  saturation  is 
expected  to  decrease  as  the  intensity  of  the  storm  increases.  To  test  this  hypothesis,  the 
same  Fort  Collins  extreme  event  used  in  the  storm  sewer  section  was  simulated  on  Dead 
Run.  The  effect  of  initial  soil  saturation  was  greatly  reduced,  as  the  saturated  test  case 
produced  a  peak  within  5%  of  the  driest  case. 
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Figure  39  Minimal  influence  of  antecedent  moisture  for  extreme  event 
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VI.  SUMMARY 


Conclusions 

Referring  to  the  six  key  objeetives  set  forth  in  the  thesis  introduetion,  the  corresponding 
conclusions  are  as  follows: 

1)  A  storm  drainage  model  was  developed  based  on  an  existing  algorithm.  It  was 
tested  for  accuracy  under  a  number  of  scenarios,  including  reverse  flow, 
backwater  effects,  and  looped  networks.  This  model  was  then  linked  to  an 
existing  distributed  hydrologic  model  to  produce  a  robust  combination  capable  of 
simulating  the  complexities  of  an  urban  watershed.  Storm  sewers  should  not  be 
modeled  as  open  channels,  since  they  allow  too  much  lateral  inflow  and  do 
represent  realistic  intake  structures  or  conveyance  properties. 

2)  Subsurface  storm  drainage  networks  have  a  significant  impact  on  flood  peaks  for 
moderate  intensity  storms.  Their  relative  importance  is  reduced,  however,  as 
rainfall  intensity  increases  and  runoff  overwhelms  the  intake  structures.  For 
moderate  storms,  storm  sewers  were  shown  to  increase  peaks  by  30%.  For 
extreme  events,  the  influence  of  storm  sewers  on  flood  magnitude  disappears. 

3)  Impervious  areas  play  an  important  role  in  the  timing  of  flood  peaks.  The 
reduction  in  surface  roughness  associated  with  parking  lots  and  roadways  reduces 
the  time  to  peak  of  moderate  storms.  However,  for  this  set  of  simulations,  the 
relative  effect  of  impervious  areas  on  flood  magnitude  is  less  than  that  of  the  total 
storm  drainage  network. 

4)  Increasing  drainage  density  has  a  strong  effect  on  flood  peaks,  but  no  influence  on 
flood  timing.  There  is  a  range  of  density  values  to  which  flood  peaks  are  very 
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sensitive.  Above  this  range,  the  network’s  effeet  on  peaks  levels  out.  Therefore, 
for  a  given  watershed  of  some  set  landuse  and  rainfall  event,  the  flood  magnitude 
approaehes  an  upper  limit. 

5)  The  width  funetion  of  a  ehannel  network  with  uniform  drainage  density  does  not 
provide  an  indieation  of  hydrologie  response.  For  non-uniform  spatial 
distribution,  the  width  funetion  was  shown  to  have  little  effeet  on  flood  peaks  for 
small  basins.  Drainage  distribution  eloser  to  the  outlet  slightly  reduees  the  time  to 
peak  diseharge.  Conveyanee  of  open  ehannel  networks  far  surpass  elosed  eonduit 
links,  and  thus  the  width  funetion  of  subsurfaee  flow  segments  should  take  into 
eonsideration  segment  geometry. 

6)  Anteeedent  soil  moisture  plays  a  very  important  role  in  flood  peaks  for  moderate 
intensity  events.  However,  its  effeet  diminishes  as  an  extreme  storm’s 
preeipitation  intensity  dwarfs  the  basin’s  infiltration  eapaeity. 
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Engineering  Recommendations 

Conclusions  from  this  research  should  not  be  blindly  used  as  guidelines  for 
development  purposes.  It  is  important  to  note  that  these  findings  pertain  to  small  basins 
of  similar  geology  and  soil  oharaoteristies.  However,  it  is  possible  to  make  some 
generalized  reeommendations  from  this  modeling  experienee.  The  dominant  eause  of 
flooding  has  been  shown  to  be  the  ability  of  a  watershed  to  quiekly  drain  its  infiltration 
exeess.  This  leads  to  two  primary  suggestions  regarding  infiltration  and  drainage 
network  properties. 

It  is  eommon  to  install  effieient  drainage  sueh  as  subsurfaee  eonerete  pipes  and 
straight,  elean  open  ehannels  to  alleviate  loealized  flooding  problems.  This  proeess 
simply  eompounds  flood  magnitudes  downstream.  Inereasing  ehannel  roughness  through 
natural  means  would  have  a  signifieant  impaet  on  attenuating  flows.  This  is  possible 
through  the  use  of  shallow,  wide,  sinuous,  grassy  swales.  Inereasing  the  storage  eapaeity 
of  these  “natural”  ehannels  eould  likewise  slow  diseharge. 

Inereasing  the  overall  quantity  of  infiltration  would  deerease  the  volumes  that  the 
ehannels  must  ultimately  handle.  The  effeet  of  impervious  areas  and  eompaeted  soils 
typieally  found  in  an  urban  setting  must  be  mitigated  by  engineered  infiltration  deviees. 
These  might  inelude  subsurfaee  storage  eavities  that  slowly  release  their  eontent  after  the 
storm,  as  well  as  eonverting  lawn  areas  to  naturally  rough  land  use  sueh  as  woodland.  It 
is  important  to  realize  that  modifieations  to  the  natural  system  have  eaused  the  inerease  in 
peak  diseharge,  and  that  re-introdueing  these  natural  meehanisms  a  way  to  reverse  the 
trend  of  urban  flooding. 
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Future  Research 


The  research  conducted  throughout  this  thesis  spawned  a  number  of  ideas  for 
future  studies.  A  number  of  these  suggestions  relate  to  the  data  of  Dead  Run  and 
Hurricane  Isabel.  It  would  be  valuable  to  model  the  watershed  with  a  finer  10  meter 
DEM  to  determine  if  the  losses  from  aggregation  were  negligible.  Improving  the  detail 
of  the  sub-surface  drainage  network  could  provide  further  analysis  of  its  influence.  As 
discussed  in  the  model  results  section,  the  single  radar  bias  value  did  not  accurately 
represent  the  second  pulse  of  rainfall.  Different  bias  values  could  be  applied  to  the  two 
storm  pulses,  while  still  conserving  storm  total  rainfall.  Similarly,  time  series  rain  gage 
data  could  be  compared  to  the  radar  results.  These  improvements  of  model  input  would 
enable  a  better  match  of  the  simulated  and  observed  outlet  hydrographs. 

Further  investigation  on  the  impact  of  impervious  areas  on  flood  timing  would 
provide  valuable  conclusions.  The  approach  might  be  to  vary  the  percent  of  impervious 
area  as  well  as  its  distribution  within  the  watershed  in  the  same  manner  as  the  drainage 
density  trials  were  executed. 

Additional  research  could  focus  on  comparing  multiple  watersheds  and  their 
drainage  densities.  By  comparing  the  upper  threshold  of  flood  magnitude  found  in  this 
thesis  to  that  of  other  basins,  one  could  determine  whether  the  drainage  density  values  are 
transferable. 
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VIII.  APPENDICES 


APPENDIX  A 

GSSHA  &  SUPERLINK  Link  Process 

The  combined  GSSHA  and  SUPERLINK  codes  will  operate  in  the  following  manner: 

1)  Initial  drainage  network  setup 

a.  GSSHA  reads  coordinates  of  all  manholes,  grates,  and  junctions  from  the 
SUPERLINK  input  fde 

b.  GSSHA  determines  the  cell  index  for  each  manhole,  grate,  and  junction 

c.  GSSHA  reads  the  channel  node/link  index  for  each  downstream  junction 
emptying  into  a  channel 

2)  GSSHA  calls  SUPERLINK 

a.  GSSHA  passes  the  water  surface  elevation  at  each  junction,  whether  in  a 
channel  or  on  the  overland  flow  plain 

b.  GSSHA  passes  the  depth  of  water  in  each  cell  containing  a  grate 

3)  SUPERLINK  is  executed 

a.  SUPERLINK  sets  boundary  conditions  for  all  downstream  junctions  with 
associated  water  surface  elevations 

b.  SUPERLINK  determines  the  depth  of  water  at  each  grate 

c.  If  the  head  at  a  grate  is  less  than  the  ground  surface  elevation, 
SUPERLINK  inserts  all  flow  into  the  grate 

d.  SUPERLINK  searches  for  any  grate,  manhole,  or  junction  heads  greater 
than  the  ground  surface  elevation  and  calculates  the  excess  to  return  to  the 
overland  plain 

e.  SUPERLINK  runs  one  timestep 

4)  SUPERLINK  returns  values  to  GSSHA 

a.  SUPERLINK  passes  the  amount  of  water  taken  from  or  added  to  an 
overland  cell 

b.  SUPERLINK  passes  the  discharge  from  downstream  junctions  to  a 
channel  node/link  or  overland  flow  cell 
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APPENDIX  B 


SUPERLINK  Input  File 

The  assimilation  of  these  three  separate  files  required  a  separate  pieee  of  eode. 
This  simple  algorithm  looped  through  eaeh  superlink,  aequiring  eonneetivity  information, 
node  positions,  and  pipe  sizes.  Quality  eontrol  was  a  serious  coneern  at  this  point,  as  the 
volume  of  manual  data  entry  left  ample  room  for  error.  Counters  were  eoded  into  the 
routine  to  seareh  for  ineorreet  number  of  oeeurrenees  of  junetions.  These  eheeks  proved 
invaluable  for  deteeting  human  blunders.  Using  the  UTM  eoordinates  of  eaeh  point,  the 
length  was  ealeulated  for  eaeh  pipe  segment.  The  output  from  this  program  was  the 
format  required  by  SUPERLINK,  as  ean  be  seen  in  this  sample.  Fields  represent 


parameters  as  defined  by  the  following  eards; 

CONNECT  "SUPERLINK  #"  "UPSTREAM  JUNCTION"  "DOWNSTREAM  JUNCTION" 

SJUNC  "JUNCTION  #"  "GROUND  SURFACE  ELEVATION  m"  "JUNCTION  BOTTOM 
ELEVATION  m"  "SURFACE  AREA  m2"  "INLET  CODE"  "UTM  NORTHING" 
"UTM  EASTING" 

SLINK  "SUPERLINK  #"  "NUMBER  OF  NODES" 

NODE  "NODE  #"  "GROUND  SURFACE  ELEVATION  m"  "INVERT  ELEVATION  m" 

"SURFACE  AREA  m2"  "INLET  CODE"  "UTM  NORTHING"  "UTM  EASTING" 

PIPE  "PIPE  #"  "SECTION  TYPE"  "DIAMETER  m"  "SLOPE"  "ROUGHNESS  n" 
"LENGTH  m" 


CONNECT 

SJUNC 

1 

1  1 

143.59 

2 

139.43 

2.000000 

1 

348491 . 615300 

4354294 . 986100 

SJUNC 

2 

137.38 

136.43 

4 . 000000 

999 

348759.887900 

4354212.379000 

SLINK 

NODE 

1 

1 

11 

143.59 

139.43 

0.000000 

777 

348491 . 615300 

4354294 . 986100 

NODE 

2 

141.73 

139.04 

1 . 000000 

2 

348528.220900 

4354281 . 109000 

NODE 

3 

143.94 

138 .81 

0.500000 

0 

348538.205500 

4354261 . 139800 

NODE 

4 

142.27 

138.46 

0.500000 

0 

348572.438500 

4354251 . 868400 

NODE 

5 

141.05 

138 . 17 

0.500000 

1 

348580.996700 

4354279.682600 

NODE 

6 

141 . 90 

137.74 

1 . 000000 

2 

348622.361500 

4354268 . 984800 

NODE 

7 

140.89 

137.48 

0.500000 

0 

348613.803300 

4354244 . 023300 

NODE 

8 

140.29 

137.03 

0.500000 

0 

348658 . 020900 

4354234.751900 

NODE 

9 

140.29 

136.88 

0.500000 

0 

348654.454900 

4354220.488200 

NODE 

10 

137.50 

136.46 

0.500000 

1 

348705.804400 

4354211.216700 

NODE 

11 

137.50 

136.45 

0.000000 

777 

348706.517600 

4354228.333200 
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NODE 

12 

137.38  136 

i.43  0.000000 

PIPE 

I 

1 

0.457200 

0.010000 

o 

o 

PIPE 

2 

1 

0.457200 

0.010000 

o 

o 

PIPE 

3 

1 

0.457200 

0.010000 

o 

o 

PIPE 

4 

1 

0.457200 

0.010000 

o 

o 

PIPE 

5 

1 

0.457200 

0.010000 

o 

o 

PIPE 

6 

1 

0.457200 

0.010000 

o 

o 

PIPE 

7 

1 

0.457200 

0.010000 

o 

o 

PIPE 

8 

1 

0.457200 

0.010000 

o 

o 

PIPE 

9 

1 

0.457200 

0.008000 

o 

o 

PIPE 

10 

1 

0.457200 

0.000500 

o 

o 

PIPE 

11 

1 

0.457200 

0.000500 

o 

o 

Ill  348759.887900  4354212.379000 

12000  39.15 

12000  22.33 

12000  35.47 

12000  29.10 

12000  42.73 

12000  26.39 

12000  45.18 

12000  14.70 

12000  52.18 

12000  17.13 

12000  55.70 
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